Vignette of analysis: Early-life demographic processes do not drive adult sex ratio biases and mating systems in sympatric coucal species

Supplementary Material A

Authors
Affiliations

Behavioural Genetics and Evolutionary Ecology, Max Planck Institute for Biological Intelligence

Department of Ornithology, Max Planck Institute for Biological Intelligence

Ignas Safari

Department Behavioural Neurobiology, Max Planck Institute for Biological Intelligence

Coucal Project

Department of Biology, College of Natural and Mathematical Sciences, University of Dodoma

Poyo Makomba

Coucal Project

Anne Hertel

Department Biology II, Ludwig Maximilians University Munich

Department Behavioural Neurobiology, Max Planck Institute for Biological Intelligence

Coucal Project

Department Biology II, Ludwig Maximilians University Munich

Published

November 24, 2023

In this document we provide all the necessary code for reproducing the analyses presented in our manuscript. To access the dataset and Rmarkdown file, please download this GitHub repository. An explanation of the files in the repository can be found in the Readme file. Please don’t hesitate to contact Luke at (luke.eberhart[at]bi.mpg.de) if you have any questions regarding the code below. For questions regarding the data collection and the study system, please contact Wolfgang (wgoymann[at]bi.mpg.de) or Safari (safariignas[at]yahoo.co.uk)

Prerequisites

GitHub repository

  • Simply follow the link to the GitHub and click on Download ZIP on the right-hand side of the page. After unzipping the folder, open the coucal_demography.Rproj file in RStudio, which will initiate a new RStudio session with the project directory referencing all the other files downloaded in the zipped folder (e.g., output/bootstraps/sensitivity_analysis/.

R packages

  • The following packages are needed for analysis and can be easily installed from CRAN or GitHub by running the following code chunk:
Code
# a vector of all the packages needed in the project
packages_required_in_project <- c("RMark",
                                  "tidyverse",
                                  "readxl",
                                  "BaSTA",
                                  "pbapply",
                                  "RColorBrewer",
                                  "grid",
                                  "Rmisc",
                                  "gss",
                                  "arm",
                                  "partR2",
                                  "parameters",
                                  "standardize",
                                  "colorBlindness",
                                  "ggthemes",
                                  "patchwork",
                                  "gt",
                                  "rptR",
                                  "tidybayes",
                                  "broom.mixed",
                                  "effects",
                                  "patchwork",
                                  "devtools",
                                  "unmarked",
                                  "R2ucare",
                                  "marked",
                                  "merTools",
                                  "bootpredictlme4",
                                  "extrafont",
                                  "survminer",
                                  "bdsmatrix",
                                  "coxme",
                                  "epitools",
                                  "survival",
                                  "magrittr",
                                  "ggpubr",
                                  "reshape2",
                                  "sp",
                                  "adehabitatLT",
                                  "reshape",
                                  "coefplot",
                                  "mapview",
                                  "lubridate",
                                  "effects",
                                  "merDeriv",
                                  "remotes")
                                  
# of the required packages, check if some need to be installed
new.packages <- 
  packages_required_in_project[!(packages_required_in_project %in% 
                                   installed.packages()[,"Package"])]

# install all packages that are not locally available
if(length(new.packages)) install.packages(new.packages)

# install the bootpredictlme4 package directly from GitHub
remotes::install_github("RemkoDuursma/bootpredictlme4")

# load all the packages into the current R session
lapply(packages_required_in_project, require, character.only = TRUE)

Plotting themes

  • The following plotting themes, colors, and typefaces are used throughout the project:
Code
# Find fonts from computer that you want. Use regular expressions to do this
# For example, load all fonts that are 'candara' or 'Candara'
extrafont::font_import(pattern = "[V/v]erdana", prompt = FALSE) 

# check which fonts were loaded
extrafont::fonts()
extrafont::fonttable()
extrafont::loadfonts() # load these into R

# define the plotting theme to be used in subsequent ggplots
luke_theme <- 
  theme_bw() +
  theme(
    text = element_text(family = "Verdana"),
    legend.title = element_text(size = 16),
    legend.text = element_text(size = 12),
    axis.title.x = element_text(size = 12),
    axis.text.x  = element_text(size = 10), 
    axis.title.y = element_text(size = 12),
    axis.text.y = element_text(size = 10),
    strip.text = element_text(size = 12),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.ticks = element_line(size = 0.5, colour = "grey40"),
    axis.ticks.length = unit(0.2, "cm"),
    panel.border = element_rect(linetype = "solid", colour = "grey"),
    legend.position = c(0.1, 0.9)
  )

# set plotting color palettes
sex_pal2 <- 
  c(pull(ggthemes_data$wsj$palettes$colors6[3,2]), 
    pull(ggthemes_data$wsj$palettes$colors6[2,2]))

sex_pal3 <- 
  c(pull(ggthemes_data$wsj$palettes$colors6[3,2]), 
    pull(ggthemes_data$wsj$palettes$colors6[3,2]),
    pull(ggthemes_data$wsj$palettes$colors6[2,2]),
    pull(ggthemes_data$wsj$palettes$colors6[2,2]))

# specify the facet labels for each species and sex
species_names <- c(
  'BC' = "Black coucal",
  'WBC' = "White-browed coucal")

sex_names <- c(
  'female' = "Females",
  'male' = "Males")

analysis_names <- c(
  'male' = "Male Mo scenario",
  'female' = "Female Mo scenario"
)

# color of mean estimate point in forest plots
col_all <- "#2E3440"

# custom color palette for the plotting of Juvenile and Adult stats
cbPalette_LTRE <- 
  c("#D9D9D9", "#D9D9D9", "#D9D9D9", 
    "#D9D9D9", "#A6A6A6", "#A6A6A6",
    "#A6A6A6")

cbPalette_sex_diff <- 
  c("#D9D9D9", "#D9D9D9", "#D9D9D9", 
    "#D9D9D9", "#A6A6A6")

# plot the comparative LTRE results
vital_rate_theme <- 
  theme_bw() +
  theme(
    text = element_text(family = "Verdana"),
    legend.position = "none",
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.ticks.length = unit(0.1, "cm"),
    panel.border = element_blank(),
    panel.spacing.x = unit(0.3, "lines"),
    panel.spacing.y = unit(0.7, "lines"),
    strip.background = element_blank()
  )
species.labs <- c("Black Coucal", "White-browed Coucal")
names(species.labs) <- c("BC", "WBC")

Analytical Functions

The following custom functions are used to bootstrap the various datasets and run the stochastic demographic model. To reproduce the results presented in our manuscript, these functions need to be stored in the environment of your R Session.

adult_CJS_bootstrap_data()

This function randomly samples the adult mark-recapture dataset to produce a random one the same size as the original, with replacement

arguments:

  • adult: the adult mark-recapture dataset
  • num_boot: a unique number for the pbapply simulation (don’t specify)
  • species: species name (either “BC” or “WBC”)
  • iter_add: incase of parallel simulations on multiple R Sessions, add a unique number here to prevent overwriting output data on disk
Code
adult_CJS_bootstrap_data <- 
  function(adult, num_boot, species, iter_add) {
  
  # sample a new adult mark-recapture dataset the same size as the original, 
  # with replacement
  adult_boot <-
    adult %>%
    sample_n(nrow(.), replace = TRUE)
  
  # make a list of the new dataset, iteration, species
  out <- list(adult_boot = adult_boot,
              iter = num_boot + ((iter_add - 1) * niter),
              species = species)
}

bootstrap_adult_survival_CJS()

runs the CJS survival analyses using the the bootstrapped sample created from adult_CJS_bootstrap_data().

arguments:

  • coucal_boot_list: the bootstrapped adult mark-recapture dataset produced from the adult_CJS_bootstrap_data() function above
  • num_boot: a unique number for the pbapply simulation (don’t specify)
  • first_year: first year of the study (i.e., the year of the first event in the encounter history)
  • bootstrap_name: a unique string for naming of files in the output directory
  • species: species name (either “BC” or “WBC”)
  • iter_add: incase of parallel simulations on multiple R Sessions, add a unique number here to prevent overwriting output data on disk
  • prefix_number: a unique prefix to prevent overwriting output data on disk (e.g., “boot_one_”)
  • out_directory: location in the project directory where the output files should be stored (e.g., “tmp”)
Code
bootstrap_adult_survival_CJS <- 
  function(coucal_boot_list, 
           num_boot, 
           first_year,
           bootstrap_name,
           species,
           iter_add,
           prefix_number,
           out_directory = "/output/bootstraps/raw/") {
    
    # specify the bootstrapped data samples (from the previous function)
    adult_ch <- coucal_boot_list[["adult_boot"]]
    
    # clean up capture histories
    adult_ch <-    
      adult_ch %>% 
      ungroup() %>% 
      as.data.frame() %>% 
      mutate(sex = as.factor(sex))
    
    # process the adult data as a "CJS"  analysis
    coucal_adult.proc <- RMark::process.data(adult_ch, 
                                             model = "CJS",
                                             groups = c("sex"),
                                             begin.time = first_year)
    
    # Create the design matrix from the processed mark-recapture datasets
    coucal_adult.ddl <- RMark::make.design.data(coucal_adult.proc)
    
    adult_survival_analysis_run = 
      function(proc_data, design_data){
        # sex- and stage-specific survival:
        Phi.sex = list(formula = ~ sex) 
        
        # Models exploring variation in encounter probability
        # constant:
        p.dot = list(formula =  ~ 1)
        
        # sex-dependent:
        p.sex = list(formula =  ~ sex)
        
        # factorial variation across year:
        p.Time = list(formula =  ~ time)
        
        # additive effects of sex and factorial year:
        p.sex_Time = list(formula =  ~ sex + time)
        
        # create a list of candidate models for all the a models above that begin with 
        # either "Phi." or "p."
        cml <-  RMark::create.model.list("CJS")
        
        # specify the data, design matrix, delete unneeded output files, and 
        # run the models in Program MARK
        model.list <-  RMark::mark.wrapper(cml, data = proc_data, 
                                           ddl = design_data, delete = TRUE, 
                                           wrap = FALSE, threads = 1, brief = TRUE,
                                           silent = TRUE, output = FALSE, prefix = prefix_number)
        
        # output the model list and sotre the results
        return(model.list)
      }
    
    # Run the models on the bootstrapped data
    adult_survival_analysis_out <-
      adult_survival_analysis_run(proc_data = coucal_adult.proc,
                                  design_data = coucal_adult.ddl)
    
    extract_top_model_output <- 
      function(rmark_output, top_model = TRUE, mod_num){
        # Find the model number for the first ranked model of the AIC table
        if(top_model == TRUE){
          mod_num <- 
            as.numeric(rownames(rmark_output$model.table[1,]))
        }
        
        else{
          mod_num <- mod_num
        }
        
        # extract and wrangle reals from model output 
        reals <- 
            rmark_output[[mod_num]]$results$real %>% 
            bind_cols(data.frame(str_split_fixed(rownames(.), " ", n = 5)), .) %>% 
            filter(X1 == "Phi") %>% 
            mutate(sex = as.factor(ifelse(unlist(str_extract_all(X2,"[FM]")) == "F", "Female", "Male"))) %>% 
            dplyr::select(sex, estimate) %>% 
            mutate(iter = num_boot + ((iter_add - 1) * niter),
                   species = species)
        
        # extract and wrangle the beta estimates from the model output
        betas <- 
          rmark_output[[mod_num]]$results$beta %>% 
          mutate(statistic = rownames(.),
                 iter = num_boot + ((iter_add - 1) * niter),
                 species = species) %>% 
          mutate(parameter = ifelse(grepl(x = statistic, pattern = "S:"), "S",
                                    ifelse(grepl(x = statistic, pattern = "p:"), "p",
                                           ifelse(grepl(x = statistic, pattern = "r:"), "r",
                                                  ifelse(grepl(x = statistic, pattern = "F:"), "F", 
                                                         ifelse(grepl(x = statistic, pattern = "Phi:"), "Phi", "XXX"))))),
                 variable = ifelse(grepl(x = statistic, pattern = "Intercept"), "Intercept",
                                   ifelse(grepl(x = statistic, pattern = "sexM:Cubic"), "sexM:Cubic",
                                          ifelse(grepl(x = statistic, pattern = "sexM:Quadratic"), "sexM:Quadratic",
                                                 ifelse(grepl(x = statistic, pattern = "sexM:Time"), "sexM:Time",
                                                        ifelse(grepl(x = statistic, pattern = "sexM"), "sexM",
                                                               ifelse(grepl(x = statistic, pattern = "Time"), "Time",
                                                                      ifelse(grepl(x = statistic, pattern = "Cubic"), "Cubic", 
                                                                             ifelse(grepl(x = statistic, pattern = "Quadratic"), "Quadratic","XXX"))))))))) %>% 
          dplyr::select(-statistic)
        
        # consolidate the AIC model selection results
        AIC_table <- 
          rmark_output$model.table %>% 
          mutate(iter = num_boot + ((iter_add - 1) * niter),
                 species = species,
                 model_no_orig = as.numeric(rownames(.))) %>% 
          mutate(model_no_rank = as.numeric(rownames(.)))
        
        # consolidate the all output into a list
        survival_model_output_list <- 
          list(reals = reals,
               betas = betas,
               AIC_table = AIC_table)
        
        survival_model_output_list
      }
    
    # extract and format survival rates from model output
    adult_survival_model_output_list <- 
      extract_top_model_output(rmark_output = adult_survival_analysis_out)
    
    # make a list of all the results from this iteration
    bootstrap_results_list <- 
      list(adult_survival_model_output = adult_survival_model_output_list,
           bootstrapped_data = coucal_boot_list)
    
    # extract file directory
    file_directory <- paste0(getwd(), "/", out_directory, "/", bootstrap_name, "_", num_boot, ".Rds")
    
    # save the results of the iteration
    save(bootstrap_results_list, file = file_directory)
    
    bootstrap_results_list
  }

run_bootstrap_adult_survival_CJS()

run the bootstrap and the CJS analysis (passed through the pbapply::pbsapply() function), see arguments defined above

arguments:

  • adult: the adult mark-recapture dataset
  • num_boot: a unique number for the pbapply() simulation (don’t specify)
  • species: species name (either “BC” or “WBC”)
  • iter_add: incase of parallel simulations on multiple R Sessions, add a unique number here to prevent overwriting output data on disk
  • first_year: first year of the study (i.e., the year of the first event in the encounter history)
  • bootstrap_name: a unique string for naming of files in the output directory
  • prefix_number: a unique prefix to prevent overwriting output data on disk (e.g., “boot_one_”)
  • out_directory: location in the project directory where the output files should be stored (e.g., “tmp”)
Code
run_bootstrap_adult_survival_CJS <- function(adult, num_boot, species, iter_add,
                                     first_year, bootstrap_name, prefix_number,
                                     out_directory)
{
  # run the sampling function and specify the datasets
  bootstrap_data_list <- 
    adult_CJS_bootstrap_data(adult = adult, 
                             num_boot = num_boot,
                             species = species, 
                             iter_add = iter_add)
  
  # run the survival analysis and ASR deduction on the sampled data
  result <- bootstrap_adult_survival_CJS(coucal_boot_list = bootstrap_data_list, 
                                         num_boot = num_boot, 
                                         first_year = first_year,
                                         bootstrap_name = bootstrap_name,
                                         species = species,
                                         iter_add = iter_add,
                                         prefix_number = prefix_number,
                                         out_directory = out_directory)
  result
}

surv_boot_out_wrangle()

loads all the output rds files from the parallel adult CJS bootstrapping procedure and consolidates them into a single object for further analysis

arguments:

  • species: species name (either “BC” or “WBC”)
  • niter: number of iterations specified in the stochastic simulation (used to predict how many files to merge from the output directory)
  • out_directory: location in the project directory where the output files from the simulation were stored (e.g., “tmp”)
Code
# boot_out_wrangle() 
surv_boot_out_wrangle <- function(species, niter, output_directory = "output/bootstraps/"){
  
  # specify directory string
  boot_out_path <- paste0(output_directory, 
                          species,"_adult_survival_CJS_bootstrap_result.rds")
  
  # load the bootstrap output file
  survival_bootstrap_result <- readRDS(boot_out_path)
  
  # bind all the adult rates output together
  adult_reals_survival_rates_boot <-
    rbind(do.call(rbind, lapply(seq(from = 1, to = niter, by = 1),
                                function(x) survival_bootstrap_result["adult_survival_model_output", x]$adult_survival_model_output$reals))) %>%
    arrange(iter, sex)
  
  # bind all the adult beta estimates output together
  adult_betas_survival_rates_boot <-
    rbind(do.call(rbind, lapply(seq(from = 1, to = niter, by = 1), 
                                function(x) survival_bootstrap_result["adult_survival_model_output", x]$adult_survival_model_output$betas)))
  
  # bind all the adult AIC table output together
  adult_AIC_tables_boot <-
    rbind(do.call(rbind, lapply(seq(from = 1, to = niter, by = 1), 
                                function(x) survival_bootstrap_result["adult_survival_model_output", x]$adult_survival_model_output$AIC_table)))
  
  # tidy all the bound data together into a mega list of results
  boot_out_tidy <- 
    list(adult_reals_survival_rates_boot = adult_reals_survival_rates_boot,
         adult_betas_survival_rates_boot = adult_betas_survival_rates_boot,
         adult_AIC_tables_boot = adult_AIC_tables_boot)
  
  boot_out_tidy
}

approx_sd()

derive the approximate sd given the upper and lower 95% confidence intervals (see:Cochrane Chapter 7.7.3.2)

arguments:

  • x1: lower bound of the confidence interval
  • x2: upper bound of the confidence interval
Code
# function to approximate the sd given the 95% CIs
approx_sd <- function(x1, x2){
  (x2-x1) / (qnorm(0.975) - qnorm(0.025))
}

matrix_ASR()

calculates the ASR of the population based on the two-sex two-stage projection matrix built by the coucal_matrix() function (below).

arguments:

  • M: a two-sex x by x projection matrix produced by the coucal_matrix() function above
  • n: x-lengthed vector representing starting stage distribution (the default is a vector with 10 individuals in each stage)
  • h: harem size index (default is 1, monogamy; polyandry < 1, polygyny > 1)
  • k: clutch size (default is 4 eggs, mean for both coucal species)
  • iterations: number of time steps to project the demogrphic model (n.b., after 100 steps our model reaches the stable stage distribution)
  • HSR: hatching sex ratio (expressed as the proportion of male chicks)
  • num_boot: a unique number for the pbapply simulation (don’t specify)
  • species: species name (either “BC” or “WBC”)
  • iter_add: incase of parallel simulations on multiple R Sessions, add a unique number here to prevent overwriting output data on disk (if you run this on one session, then don’t specify)
Code
matrix_ASR <-
  function(M, 
           n = rep(10, nrow(M)), 
           h = 1, 
           k = 4, 
           iterations = 1000, 
           HSR = 0.5, 
           num_boot, 
           species, 
           iter_add){
    
    # Number of stages in matrix
    x <- length(n)
    
    # Number of time steps to simulate
    t <- iterations
    
    # an empty t by x matrix to store the stage distributions
    stage <- matrix(data = numeric(x * t), nrow = x, ncol = t)
    
    # an empty t by x matrix to store the stage-specific frequencies
    pop_dist <- matrix(data = numeric(x * t), nrow = x, ncol = t)
    rownames(pop_dist) <- rownames(M)
    colnames(pop_dist) <- 0:(t - 1)
    
    # an empty t vector to store the population sizes
    pop <- numeric(t)
    
    # for loop that goes through each of t time steps
    for (i in 1:t) {
      
      # stage distribution at time t
      stage[,i] <- n
      
      # population size at time t
      pop[i] <- sum(n)
      
      # number of male adults at time t = (number alive at t-1) * (survival rate)
      M2 <- stage[4, i] * M["M_Adult", "M_Adult"]
      
      # number of female adults at time t
      F2 <- stage[2, i] * M["F_Adult", "F_Adult"]
      
      # Female freq-dep fecundity of Female chicks
      M["F_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h)) * (1 - HSR))
      
      # Female freq-dep fecundity of Male chicks
      M["M_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h)) * HSR)
      
      # Male freq-dep fecundity of Female chicks
      M["F_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h)) * (1 - HSR))
      
      # Male freq-dep fecundity of Male chicks
      M["M_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h)) * HSR)
      
      # define the new n (i.e., new stage distribution at time t)
      n <- M %*% n
      
      # define rownames of stage matrix
      rownames(stage) <- rownames(M)
      
      # define colnames of stage matrix
      colnames(stage) <- 0:(t - 1)
      
      # calculate the proportional stable stage distribution
      stage <- apply(stage, 2, function(x) x/sum(x))
      
      # define stable stage as the last stage
      stable.stage <- stage[, t]
      
      pop_dist[, i] <- n
      
    }
    # calc ASR as the proportion of the adult stable stage class that is male
    ASR <- stable.stage[4]/(stable.stage[2] + stable.stage[4])
    SSR <- stable.stage[3]/(stable.stage[1] + stable.stage[3])
    
    # make a list of results
    pop.proj <- list(ASR = ASR,
                     SSR = SSR,
                     lambda = pop[t]/pop[t - 1],
                     pop_size = pop,
                     stage_size = pop_dist,
                     stable.stage = stable.stage,
                     stage.vectors = stage,
                     SSD_M2 = stable.stage[4],
                     SSD_F2 = stable.stage[2],
                     iter = num_boot + ((iter_add - 1) * niter), 
                     species = species)
    
    # print the list as output to the function
    pop.proj
  }

coucal_matrix()

builds the two-sex Lefkovitch matrix using the vital rates specified in the demographic_rates object. This function also give the option to specify a one-sex matrix for analytical comparisons between 2- and 1-sex matrix models

arguments:

  • demographic_rates_: a list of all the vital rates of an iteration of the stochastic simulation run within the run_bootstrap_juv_hazd_ad_surv_ASR() function below
  • two_sex: whether to make a one-sex or two-sex matrix of the rates provided (default is ‘two_sex = TRUE’)
Code
coucal_matrix <- 
  function(demographic_rates_, two_sex = TRUE){
    if(two_sex){
      
      # Define coucal life-stages of the coucal matrix model
      stages <- c("F_1st_year",  "F_Adult",  "M_1st_year",  "M_Adult")
      
      # Build the 4x4 matrix
      result <- 
        matrix(c(
          
          # top row of matrix
          0, NA, 0, NA, 
          
          # second row of matrix
          (demographic_rates_$Egg_survival * 
             demographic_rates_$F_Nestling_survival * 
             demographic_rates_$F_Groundling_survival *
             demographic_rates_$F_Fledgling_survival),
          demographic_rates_$F_Adult_survival,
          0, 0,
          
          # third row of matrix
          0, NA, 0, NA, 
          
          # fourth row of matrix
          0, 0, 
          (demographic_rates_$Egg_survival * 
             demographic_rates_$M_Nestling_survival * 
             demographic_rates_$M_Groundling_survival *
             demographic_rates_$M_Fledgling_survival),
          demographic_rates_$M_Adult_survival),
          nrow = length(stages), byrow = TRUE,
          dimnames = list(stages, stages))
    }
    else{
      
      # Define coucal life-stages of the Ceuta snowy coucal matrix model
      stages <- c("1st_year",  "Adult")
      
      # Build the 4x4 matrix
      result <- 
        matrix(c(
          
          # top row of matrix
          0, NA,
          
          # second row of matrix
          (demographic_rates_$Egg_survival * 
             demographic_rates_$Nestling_survival * 
             demographic_rates_$Groundling_survival *
             demographic_rates_$Fledgling_survival),
          demographic_rates_$Adult_survival),
          nrow = length(stages), byrow = TRUE,
          dimnames = list(stages, stages))
    }
    result
  }

bootstrap_hazard_data()

This function randomly samples one sibling per nest from the offspring survival dataset

arguments:

  • offspring: the radio-tracked offspring dataset
  • num_boot: a unique number for the pbapply simulation (don’t specify)
  • species: species name (either “BC” or “WBC”)
  • iter_add: incase of parallel simulations on multiple R Sessions, add a unique number here to prevent overwriting output data on disk
  • alpha_value: Parameter defining cross-validation score for smoothing parameter selection in the gss::sshzd() function (n.b., Gu 2014 J. Stat Software 58:art5 recommends alpha = 1.4)
  • max_time: maximum number of days since hatch to model survival (‘70’ in the case of our analysis)
  • niter: number of iterations specified in the stochastic simulation
Code
# bootstrap_hazard_data() 
bootstrap_hazard_data <- 
  function(offspring, num_boot, species, iter_add, alpha_value, max_time, niter) {
    
    # set attempt to 0 at start of each loop
    attempt <- 0
    
    max_time <- max_time
    
    hzd_ss_try_M_lambda <- 1.1
    hzd_ss_try_F_lambda <- 1.1

    # store simulated estimates only if lambda is less than 1.1
    while( (hzd_ss_try_M_lambda > -1.1 | hzd_ss_try_F_lambda > 1.1) && attempt <= 100 ) {

      # next attempt
      attempt <- attempt + 1

      # sample a new offspring dataset containing only one nest member per draw
      offspring_boot <-
        offspring %>%
        group_by(nest_ID) %>%
        sample_n(1)

      try(
        hzd_ss_try_M <- sshzd(Surv(exit, event, entry) ~ exit,
                              data = filter(offspring_boot, sex == "M"),
                              alpha = alpha_value),
        silent = TRUE
      )

      try(
        hzd_ss_try_F <- sshzd(Surv(exit, event, entry) ~ exit,
                              data = filter(offspring_boot, sex == "F"),
                              alpha = alpha_value),
        silent = TRUE
      )

      # store calculated peak (i.e., the apex of the polynomial curve)
      try(
        hzd_ss_try_M_lambda <- hzd_ss_try_M$lambda
      )
      # store calculated peak (i.e., the apex of the polynomial curve)
      try(
        hzd_ss_try_F_lambda <- hzd_ss_try_F$lambda
      )
    }
    
    # store simulated estimates only if the hzdcurve can be estimated
    while( (exists("hzd_curve_try_M") == FALSE | exists("hzd_curve_try_F") == FALSE) && attempt <= max_time) {
      
      time_vector <- seq(0, max_time - attempt, 1)
      
      # next attempt
      attempt <- attempt + 1
      
      try(
        hzd_ss_try_M <- sshzd(Surv(exit, event, entry) ~ exit, 
                              data = filter(offspring_boot, sex == "M"), 
                              alpha = alpha_value),
        silent = TRUE
      )
      
      try(
        hzd_ss_try_F <- sshzd(Surv(exit, event, entry) ~ exit, 
                              data = filter(offspring_boot, sex == "F"), 
                              alpha = alpha_value),
        silent = TRUE
      )
      
      # simulate an estimate
      try(
        hzd_curve_try_M <- 
          hzdcurve.sshzd(object = hzd_ss_try_M, time = time_vector, se = TRUE),
        silent = TRUE
      )
      try(
        hzd_curve_try_F <- 
          hzdcurve.sshzd(object = hzd_ss_try_F, time = time_vector, se = TRUE),
        silent = TRUE
      )
      
    }
    
    # make a list including elements that will be used in the next function
    out <- list(offspring_boot = offspring_boot, 
                iter = num_boot + ((iter_add - 1) * niter),
                species = species,
                time_vector = time_vector,
                lambdaM = hzd_ss_try_M_lambda,
                lambdaF = hzd_ss_try_F_lambda)
  }

run_bootstrap_juv_hazd_ad_surv_ASR()

initiates the bootstrap_hazard_data() and runs the juvenile survival analysis to derive sex- and stage-specific survival rates. Then constructs the matrix model with the rates and projects this into the future to acquire the stable stage distribution and derive the ASR.

arguments:

  • num_boot: a unique number for the pbapply simulation (don’t specify)
  • offspring: the radio-tracked offspring dataset
  • k_dist: a two-number vector containing the mean and sd of clutch size (output from the analyses below)
  • HSR_dist: a two-number vector containing the mean and sd of hatching sex ratio (output from the analyses below)
  • h_dist: a two-number vector containing the mean and sd of the harem size index (output from the analyses below)
  • egg_surv_dist: a two-number vector containing the mean and sd of egg survival (output from the analyses below)
  • flight_age_distF: a two-number vector containing the mean and sd of female-specific flight age (output from the analyses below)
  • flight_age_distM: a two-number vector containing the mean and sd of male-specific flight age (output from the analyses below)
  • fledge_age_distF: a two-number vector containing the mean and sd of female-specific fledge age (output from the analyses below)
  • fledge_age_distM: a two-number vector containing the mean and sd of male-specific fledge age (output from the analyses below)
  • bootstrap_name: a unique string for naming of files in the output directory
  • adult_surival_boot_out: the bootstrapped adult encounter history data
  • species: species name (either “BC” or “WBC”)
  • iter_add: incase of parallel simulations on multiple R Sessions, add a unique number here to prevent overwriting output data on disk
  • prefix_number: a unique prefix to prevent overwriting output data on disk (e.g., “boot_one_”)
  • alpha_value: Parameter defining cross-validation score for smoothing parameter selection in the gss::sshzd function (default is 1.4).
  • max_time: maximum number of days since hatch to model survival (‘70’ in the case of our analysis)
  • niter: number of iterations specified in the stochastic simulation
  • out.directory: location in the project directory where the output files should be stored (e.g., “tmp”)
Code
run_bootstrap_juv_hazd_ad_surv_ASR <- function(num_boot, 
                                               offspring,
                                               k_dist, 
                                               HSR_dist, 
                                               h_dist, 
                                               egg_surv_dist, 
                                               flight_age_distF,
                                               flight_age_distM,
                                               fledge_age_distF, 
                                               fledge_age_distM,
                                               bootstrap_name, 
                                               adult_surival_boot_out,
                                               species, 
                                               iter_add, 
                                               prefix_number,
                                               alpha_value = 1.4,
                                               max_time = 70,
                                               niter,
                                               output.directory)
{
  # run the sampling function and specify the datasets
  coucal_boot_list <- 
    bootstrap_hazard_data(offspring = offspring, 
                          num_boot = num_boot,
                          species = species, 
                          iter_add = iter_add,
                          alpha_value = alpha_value,
                          max_time = max_time,
                          niter = niter)
  
  # run the survival analysis and ASR deduction on the sampled data
  # specify the bootstrapped data samples (from the previous function)
    offspring_data <- coucal_boot_list[["offspring_boot"]]
    
    # clean up capture histories
    offspring_data <-    
      offspring_data %>% 
      ungroup() %>% 
      as.data.frame()
    
    # fit smoothed spline of hazard function for either sex
    M_haz_ss <- sshzd(Surv(exit, event, entry) ~ exit, 
                      data = filter(offspring_data, sex == "M"), 
                      alpha = alpha_value)

    F_haz_ss <- sshzd(Surv(exit, event, entry) ~ exit, 
                      data = filter(offspring_data, sex == "F"), 
                      alpha = alpha_value)

    haz_ss_lambda <- list(Male_haz_ss_lambda = M_haz_ss$lambda,
                          Female_haz_ss_lambda = F_haz_ss$lambda)
    
    # extract fitted estimates from the spline function
    M_haz_ss_curve <- 
      hzdcurve.sshzd(object = M_haz_ss, time = coucal_boot_list[["time_vector"]], se = TRUE)
    
    F_haz_ss_curve <- 
      hzdcurve.sshzd(object = F_haz_ss, time = coucal_boot_list[["time_vector"]], se = TRUE)
    
    haz_ss_curve <- 
      expand.grid(species = species, 
                  age = coucal_boot_list[["time_vector"]],
                  sex = c("Male", "Female")) %>% 
      mutate(fit = c(M_haz_ss_curve$fit, F_haz_ss_curve$fit),
             se = c(M_haz_ss_curve$se, F_haz_ss_curve$se)) %>% 
      mutate(estimate = 1 - fit,
             upper = 1 - fit * exp(1.96 * se),
             lower = 1 - fit / exp(1.96 * se),
             iter = num_boot)
    
    # transform the daily nestling survival (DCS) to apparent fledgling success
    # by calculating the product of all DCS estimates:
    
    fledge_ageM <- rnorm(1, fledge_age_distM[1], fledge_age_distM[2])
    fledge_ageF <- rnorm(1, fledge_age_distF[1], fledge_age_distF[2])
    
    flight_ageM <- rnorm(1, flight_age_distM[1], flight_age_distM[2])
    flight_ageF <- rnorm(1, flight_age_distF[1], flight_age_distF[2])
    
    coucal_nestling_survivalF <-
      haz_ss_curve %>% 
      filter(age <= fledge_ageF & sex == "Female") %>% 
      group_by(sex) %>%
      dplyr::summarise(value = prod(estimate), .groups = 'drop') %>%
      mutate(stage = "nestling",
             rate = "survival")
    
    coucal_nestling_survivalM <-
      haz_ss_curve %>% 
      filter(age <= fledge_ageM & sex == "Male") %>% 
      group_by(sex) %>%
      dplyr::summarise(value = prod(estimate), .groups = 'drop') %>%
      mutate(stage = "nestling",
             rate = "survival")
    
    coucal_nestling_survival <- 
      bind_rows(coucal_nestling_survivalF, coucal_nestling_survivalM)
    
    coucal_groundling_survivalF <-
      haz_ss_curve %>% 
      filter(age < (flight_ageF + fledge_ageF) & age > fledge_ageF & sex == "Female") %>% 
      group_by(sex) %>% 
      dplyr::summarise(value = prod(estimate), .groups = 'drop') %>% 
      mutate(stage = "groundling",
             rate = "survival")
    
    coucal_groundling_survivalM <-
      haz_ss_curve %>% 
      filter(age < (flight_ageM + fledge_ageM) & age > fledge_ageM & sex == "Male") %>% 
      group_by(sex) %>% 
      dplyr::summarise(value = prod(estimate), .groups = 'drop') %>% 
      mutate(stage = "groundling",
             rate = "survival")
    
    coucal_groundling_survival <- 
      bind_rows(coucal_groundling_survivalF, coucal_groundling_survivalM)
    
    coucal_fledgling_survivalF <-
      haz_ss_curve %>% 
      filter(age >= (flight_ageF + fledge_ageF) & sex == "Female") %>% 
      group_by(sex) %>% 
      dplyr::summarise(value = prod(estimate), .groups = 'drop') %>% 
      mutate(stage = "fledgling",
             rate = "survival")
    
    coucal_fledgling_survivalM <-
      haz_ss_curve %>% 
      filter(age >= (flight_ageM + fledge_ageM) & sex == "Male") %>% 
      group_by(sex) %>% 
      dplyr::summarise(value = prod(estimate), .groups = 'drop') %>% 
      mutate(stage = "fledgling",
             rate = "survival")
    
    coucal_fledgling_survival <- 
      bind_rows(coucal_fledgling_survivalF, coucal_fledgling_survivalM)
    
    coucal_adult_survival <-
      filter(adult_surival_boot_out, iter == sample(1:niter, 1)) %>%
      dplyr::select(-species, -iter) %>%
      mutate(stage = "adult",
             rate = "survival")
    
    egg_survival <- rnorm(1, egg_surv_dist[1], egg_surv_dist[2])
    
    coucal_egg_survival <- 
      data.frame(sex = NA,
                 value = egg_survival,
                 stage = c("egg"),
                 rate = c("survival"))
    
    coucal_mating_system <- 
      data.frame(sex = NA,
                 value = 1/rnorm(1, h_dist[1], h_dist[2]),
                 stage = c("h"),
                 rate = c("fecundity"))
    
    coucal_clutch_size <- 
      data.frame(sex = NA,
                 value = rnorm(1, k_dist[1], k_dist[2]),
                 stage = c("k"),
                 rate = c("fecundity"))
    
    coucal_HSR <- 
      data.frame(sex = NA,
                 value = rnorm(1, HSR_dist[1], HSR_dist[2]),
                 stage = c("HSR"),
                 rate = c("fecundity"))
    
    coucal_flight_age <- 
      data.frame(sex = c("Female", "Male"),
                 value = c(flight_ageF, flight_ageM),
                 stage = c("flight_age"),
                 rate = c("development"))
    
    coucal_fledge_age <- 
      data.frame(sex = c("Female", "Male"),
                 value = c(fledge_ageF, fledge_ageM),
                 stage = c("fledge_age"),
                 rate = c("development"))
    
    # Bind the juvenile and adult dataframe with the nestlings
    coucal_vital_rates <- 
      bind_rows(coucal_egg_survival,
                coucal_nestling_survival,
                coucal_groundling_survival,
                coucal_fledgling_survival,
                coucal_adult_survival,
                coucal_mating_system,
                coucal_clutch_size,
                coucal_HSR,
                coucal_flight_age,
                coucal_fledge_age) %>% 
      as.data.frame() %>% 
      mutate(iter = num_boot, #+ ((iter_add - 1) * niter),
             species = species) %>% 
      arrange(sex, rate, stage)
    
    # Create a list of demographic rates from the survival analyses above
    demographic_rates <- list(Egg_survival = pull(filter(coucal_vital_rates, stage == "egg"), value),
                              F_Nestling_survival = pull(filter(coucal_vital_rates, 
                                                                stage == "nestling" & rate == "survival" & sex == "Female"), 
                                                         value),
                              F_Groundling_survival = pull(filter(coucal_vital_rates, 
                                                                  stage == "groundling" & rate == "survival" & sex == "Female"), 
                                                           value),
                              F_Fledgling_survival = pull(filter(coucal_vital_rates, 
                                                                 stage == "fledgling" & rate == "survival" & sex == "Female"), 
                                                          value),
                              F_Adult_survival = pull(filter(coucal_vital_rates, 
                                                             stage == "adult" & rate == "survival" & sex == "Female"), 
                                                      value),
                              M_Nestling_survival = pull(filter(coucal_vital_rates, 
                                                                stage == "nestling" & rate == "survival" & sex == "Male"), 
                                                         value),
                              M_Groundling_survival = pull(filter(coucal_vital_rates, 
                                                                  stage == "groundling" & rate == "survival" & sex == "Male"), 
                                                           value),
                              M_Fledgling_survival = pull(filter(coucal_vital_rates, 
                                                                 stage == "fledgling" & rate == "survival" & sex == "Male"), 
                                                          value),
                              M_Adult_survival = pull(filter(coucal_vital_rates, 
                                                             stage == "adult" & rate == "survival" & sex == "Male"), 
                                                      value),
            
                              # Define hatching sex ratio
                              HSR = pull(filter(coucal_vital_rates, stage == "HSR"), value),
                              
                              # Define the mating system (h), and clutch size (k)
                              h = pull(filter(coucal_vital_rates, stage == "h"), value),
                              k = pull(filter(coucal_vital_rates, stage == "k"), value))
    
    # Build matrix based on rates specified in the list above
    demographic_matrix <- coucal_matrix(demographic_rates_ = demographic_rates)
    
    # Determine the ASR at the stable stage distribution
    ASR_SSD <- matrix_ASR(M = demographic_matrix,
                          h = demographic_rates$h,
                          HSR = demographic_rates$HSR, 
                          iterations = 100,
                          num_boot = num_boot,
                          species = species,
                          iter_add = 1,
                          n = rep(10, nrow(demographic_matrix)))
    
    # Extract ASR
    ASR_estimate <- ASR_SSD$ASR
    
    # make a list of all the results from this iteration
    bootstrap_results_list <- 
      list(offspring_hazard_lambda = haz_ss_lambda, 
           offspring_hazard_rates = haz_ss_curve,
           coucal_vital_rates = coucal_vital_rates,
           ASR_SSD = ASR_SSD,
           bootstrapped_data = coucal_boot_list)
    
    # extract file directory
    file_directory <- paste0(getwd(), output.directory, "/", bootstrap_name, "_", num_boot, ".Rds")
    
    # save the results of the iteration
    save(bootstrap_results_list, file = file_directory)
    
    result = bootstrap_results_list
  
  result
}

hazard_boot_out_wrangle()

If the product of run_bootstrap_juv_hazd_ad_surv_ASR() is saved on disk, this function will tidy up the components of the raw output for subsequent analyses

arguments:

  • species: species name (either “BC” or “WBC”)
  • niter: number of iterations specified in the stochastic simulation
  • output_dir: location in the project directory where the output files are to be found by the function (e.g., “tmp”)
  • rds_file: name of the stochastic simulation rds file saved to disk
Code
# single_model_boot_out_wrangle() 
hazard_boot_out_wrangle <- 
  function(species, niter, output_dir, rds_file){
    
    # specify directory string
    boot_path <- paste0(output_dir, species, rds_file, ".rds")
    
    # load each of the four bootstrap output files
    bootstrap_out <- readRDS(file = boot_path)
    
    # bind all the fledgling daily survival rates output together
    hazard_rates_boot <-
      do.call(rbind, lapply(seq(from = 1, to = niter, by = 1),
                            function(x) bootstrap_out["offspring_hazard_rates", x]$offspring_hazard_rates)) %>%
      arrange(iter, sex, age)
    
    # bind all the ASR estimates output together
    ASR_ests_boot <-
      do.call(rbind, lapply(seq(from = 1, to = niter, by = 1),
                            function(x) bootstrap_out["ASR_SSD", x]$ASR_SSD$ASR)) %>% 
      as.data.frame() %>% 
      mutate(species = species)
    
    # bind all the ASR estimates output together
    vital_rate_ests_boot <-
      do.call(rbind, lapply(seq(from = 1, to = niter, by = 1),
                            function(x) bootstrap_out["coucal_vital_rates", x]$coucal_vital_rates))
    
    # tidy all the bound data together into a mega list of results
    boot_out_tidy <- 
      list(hazard_rates_boot = hazard_rates_boot,
           vital_rate_ests_boot = vital_rate_ests_boot,
           ASR_ests_boot = ASR_ests_boot)
    
    boot_out_tidy
  }

sex_diff_hazard()

takes the list of results from the bootstrap procedure and calculates the sex difference in survival at each life stage for each iteration

arguments:

  • boot_out_list: the object produced by the hazard_boot_out_wrangle() above
  • niter: number of iterations specified in the stochastic simulation
Code
sex_diff_hazard <- function(boot_out_list, niter) {
  
  # make an empty datarame to store the results
  sex_diff_surv_output <- data.frame(Adult = niter,
                                     Fledgling = niter,
                                     Groundling = niter,
                                     Nestling = niter,
                                     ISR = niter,
                                     HSR = niter,
                                     h = niter)
  
  # for loop to go through each iteration and calculate the difference between 
  # female and male survival rates for each stage.
  for(i in 1:niter){
    Adult <- 
      pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], 
                  stage == "adult" & rate == "survival" & sex == "Female"), value) -
      pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], 
                  stage == "adult" & rate == "survival" & sex == "Male"), value)
    Fledgling <- 
      pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], 
                  stage == "fledgling" & rate == "survival" & sex == "Female"), value) -
      pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], 
                  stage == "fledgling" & rate == "survival" & sex == "Male"), value)
    Groundling <- 
      pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], 
                  stage == "groundling" & rate == "survival" & sex == "Female"), value) -
      pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], 
                  stage == "groundling" & rate == "survival" & sex == "Male"), value)
    Nestling <- 
      pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], 
                  stage == "nestling" & rate == "survival" & sex == "Female"), value) -
      pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], 
                  stage == "nestling" & rate == "survival" & sex == "Male"), value)
    HSR <- 
      0.5 - pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], 
                  stage == "HSR"), value)
    h <- 
      1 - pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], 
                  stage == "h"), value)
    
    sex_diff_surv_output[i, 1] <- ifelse(length(Adult) > 0, Adult, NA)
    sex_diff_surv_output[i, 2] <- ifelse(length(Fledgling) > 0, Fledgling, NA)
    sex_diff_surv_output[i, 3] <- ifelse(length(Groundling) > 0, Groundling, NA)
    sex_diff_surv_output[i, 4] <- ifelse(length(Nestling) > 0, Nestling, NA)
    sex_diff_surv_output[i, 6] <- ifelse(length(HSR) > 0, HSR, NA)
    sex_diff_surv_output[i, 7] <- ifelse(length(h) > 0, h, NA)
    
    
  }
  
  # restructure the output and lable columns
  sex_diff_surv_output <- reshape2::melt(data = sex_diff_surv_output)
  colnames(sex_diff_surv_output) <- c("stage", "difference")
  
  # return the output
  sex_diff_surv_output
}

Wrangling survival data

The following chunks show how we wrangled the raw field data to consolidate the data into the appropriate format for the survival analyses

Offspring encounter and census history

Code
status_dat_all <- 
  # read raw data
  read.csv("data/raw/Coucal_chick_survival_2001-2019_20200129.csv", 
           header = TRUE, stringsAsFactors = FALSE, na.strings = c("", " ", "NA")) %>% 
  
  # rename ring_ID column
  dplyr::rename(ring_ID = Ring_ID) %>% 
  
  # make all entries lower case for consistency
  mutate(Fledged_status = tolower(Fledged.),
         site = tolower(site)) %>% 
  
  # select variables of interest
  dplyr::select(species, ring_ID, lab_no, sex, year, site, nest_ID, pref_age, 
                Fledged_status, postf_age, postf_status, ageC, statusC, 
                lay_date, hatch_order) %>% 
  
  # remove all white space from data
  mutate(across(everything(), ~str_trim(.x))) %>% 
  mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>% 
  
  # specify empty data as NA
  mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>% 
  
  # exclude all individuals that died in the nest
  # filter(Fledged_status == "yes") %>% 
  
  # classify columns
  mutate(sex = as.factor(sex),
         ageC = as.numeric(ageC),
         postf_age = as.numeric(postf_age),
         postf_status = as.numeric(postf_status),
         hatch_order = as.numeric(hatch_order),
         pref_age = as.numeric(pref_age)) %>% 
  
  # remove rows with missing sex, age, and status info
  filter(!is.na(sex) & !is.na(ageC) & !is.na(statusC)) %>% 
  
  # make a unique id for each individual
  mutate(ind_ID = paste(nest_ID, lab_no, ring_ID, sep = "_"),
         
         # create the age of entry into the data (all at age 15)
         entry = 0,
         
         # specify the age of death or censoring
         exit = ageC,
         
         # make the event numeric and specify if 
         # the individual died (1) or was censored (0)
         event = as.numeric(statusC)) %>% 
  
  # consolidate to variables of interest
  dplyr::select(species, ind_ID, nest_ID, year, sex, entry, exit, event, hatch_order)

Adult encounter history

Code
detect_dat <-
  read_xlsx("data/raw/All_coucal_waypoints_2001_2019_20200202.xlsx", na = "NA", col_types = "text") %>% 
  dplyr::select(species, ring_ID, sex, year, site, age_status, date_dec) %>% 
  mutate(month = str_sub(date_dec, start = 5, end = 6),
         day = str_sub(date_dec, start = 7, end = 8)) %>% 
  mutate(date = as.Date(paste(year, month, day, sep = "-"), format = "%Y-%m-%d")) %>% 
  dplyr::select(-month, -day, -date_dec) %>% 
  mutate(across(c("sex", "site", "age_status"), tolower)) %>%
  mutate(sex = ifelse(sex == "female", "F", ifelse(sex == "male", "M", "XXX")),
         age_status = ifelse(age_status == "adult", "A", ifelse(age_status == "juvenile", "J", "XXX")),
         ring_ID = str_replace_all(string = ring_ID, fixed(" "), "")) %>%
  mutate(across(c("species", "ring_ID", "sex", "site", "age_status"), as.factor))

# make annual capture histories for adults
BC_detect_dat_A <-
  detect_dat %>% 
  filter(age_status == "A" & species == "BC")

# use the BaSTA function "CensusToCaptHist()" function to convert long format
# encounter histories of each individual, to wide format with 1's and 0's for 
# each year of encounter
BC_detect_dat_A_ch <- 
  CensusToCaptHist(ID = BC_detect_dat_A$ring_ID,
                   d = BC_detect_dat_A$year) %>% 
  dplyr::rename(ring_ID = ID) %>%
  mutate(ID = rownames(.)) %>% 
  left_join(., dplyr::select(BC_detect_dat_A, ring_ID, sex, year), by = "ring_ID") %>% 
  distinct()

Black_Coucal_adult_CJS_ch <- 
  data.frame(ch = apply(BC_detect_dat_A_ch[, 2:19 ] , 1, paste, collapse = "")) %>% 
  bind_cols(., dplyr::select(BC_detect_dat_A_ch, sex, ring_ID)) %>% 
  mutate(across(everything(), ~str_trim(.x))) %>% 
  mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>% 
  mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>% 
  distinct()

WBC_detect_dat_A <-
  detect_dat %>% 
  filter(age_status == "A" & species == "WBC") %>% 
  mutate(year = as.numeric(year))

WBC_detect_dat_A_ch <- 
  CensusToCaptHist(ID = WBC_detect_dat_A$ring_ID,
                   d = WBC_detect_dat_A$year) %>% 
  dplyr::rename(ring_ID = ID) %>%
  mutate(ID = rownames(.)) %>% 
  left_join(., dplyr::select(WBC_detect_dat_A, ring_ID, sex, year), by = "ring_ID") %>% 
  distinct()

White_browed_Coucal_adult_CJS_ch <- 
  data.frame(ch = apply(WBC_detect_dat_A_ch[, 2:16 ] , 1, paste, collapse = "")) %>% 
  bind_cols(., dplyr::select(WBC_detect_dat_A_ch, sex, ring_ID)) %>% 
  mutate(across(everything(), ~str_trim(.x))) %>% 
  mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>% 
  mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>% 
  distinct()

Bootstrapped mark-recapture analysis of White-browed Coucal adults

Here we resample the annual encounter history of adult White-browed Coucals and conduct a CJS mark-recapture analysis to derive the mean estimate and variance in annual apparent survival.

Set the ‘niter’ to run the bootstrap over several iterations (e.g., 10000). For reproducibility, set the seed with the set.seed() function. Set the ‘output.directory’ so that the subsequent chunks can refer to the appropriate folder within the project directory.

Code
niter = 10000
set.seed(14)
output.directory = "output/bootstraps/raw/"

run bootstrap procedure on White-browed Coucals

Code
WBC_adult_survival_CJS_bootstrap_result <-
  pbsapply(1:niter, run_bootstrap_adult_survival_CJS, 
           adult = White_browed_Coucal_adult_CJS_ch,
           species = "WBC",
           iter_add = 1,
           first_year = 2005,
           bootstrap_name = "WBC_boot_one",
           prefix_number = "boot_CJS_",
           out_directory = "output/bootstraps/raw")
Code
# save the output
saveRDS(WBC_adult_survival_CJS_bootstrap_result,
        file = "output/bootstraps/raw/WBC_adult_survival_CJS_bootstrap_result.rds")

# tidy it up
WBC_adult_survival_CJS_bootstrap_tidy <- 
  surv_boot_out_wrangle(species = "WBC", niter = niter, output_directory = output.directory)

# save the tidy version
saveRDS(WBC_adult_survival_CJS_bootstrap_tidy,
        file = "output/bootstraps/tidy/WBC_adult_survival_CJS_bootstrap_tidy.rds")

Derive the Black Coucal sex-specific survival rates by transforming White-browed Coucal sex-specific survival rates to reflect the ~0.7 ASR that is observed annually

Code
# bootstrapped WBC survival rates
WBC_adult_surival_boot_out <-
  WBC_boot_out$adult_reals_survival_rates_boot %>% 
  dplyr::rename(value = estimate)

# boostrapped BC survival rates (i.e., derived from the WBC rates)
trans_WBC_adult_surival_boot_out <-
  WBC_boot_out$adult_reals_survival_rates_boot %>% 
  dplyr::rename(value = estimate) %>% 
  mutate(value = ifelse(sex == "Female", value * 0.6, value * 1.4))

Parameterizing demographic rates

Each of the following rates has a mean and standard deviation that are used in the stochastic simulation of adult sex ratio.

Mating system

Code
mating_dat <- 
  # read raw data
  read.csv("data/raw/Coucal_Nr_mates_2001_2019_20200123.csv", 
           header = TRUE, stringsAsFactors = FALSE, na.strings = c("", " ", "NA")) %>% 

  # make all entries lower case for consistency
  mutate(site = tolower(site),
         age_status = tolower(age_status),
         sex = tolower(sex)) %>% 

  # specify empty data as NA
  mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>% 
  
  # classify columns
  mutate(sex = as.factor(sex),
         species = as.factor(species),
         Nr_partners = as.numeric(Nr_partners),
         age_status = as.factor(age_status),
         site = as.factor(site)) %>% 
  
  # remove rows with missing sex and age, and post-fledged status info
  filter(!is.na(sex) & !is.na(Nr_partners) & !is.na(species) & species != "CTC") %>% 
  
  # make a unique id for each individual
  mutate(ind_ID = str_replace_all(ring_ID, fixed(" "), "")) %>% 
  mutate(ID_partner.s. = str_replace_all(ID_partner.s., fixed(" "), "")) %>% 

  # make a unique id for each individual
  mutate(ind_ID = str_replace_all(ind_ID, fixed("_"), "")) %>% 
  mutate(ID_partner.s. = str_replace_all(ID_partner.s., fixed("_"), "")) %>% 
  
  
  # consolidate to variables of interest
  dplyr::select(year, species, ind_ID, site, sex, age_status, Nr_partners, ID_partner.s.)

sex_specific_mating_system <- 
  mating_dat %>% 
  group_by(ind_ID, sex, species) %>% 
  dplyr::summarise(mean_annual_no_mates_ind = mean(Nr_partners, na.rm = TRUE),
                   n_years = n_distinct(year)) %>% 
  mutate(mu_i = ifelse(mean_annual_no_mates_ind <= 1, 1, mean_annual_no_mates_ind)) %>% 
  group_by(species, sex) %>% 
  dplyr::summarise(mean_annual_no_mates = mean(mu_i, na.rm = TRUE),
                   var_annual_no_mates = var(mu_i, na.rm = TRUE),
                   median_annual_no_mates = median(mu_i, na.rm = TRUE),
                   sd_annual_no_mates = sd(mu_i, na.rm = TRUE),
                   n = n_distinct(ind_ID), .groups = "drop") %>% 
  mutate(sample_size = paste("n = ", n, sep = ""),
         species_plot = ifelse(species == "WBC", 1.9, 1.1),
         species_lab = ifelse(species == "WBC", 1, 2))

# To obtain a female-based h index for each population the inverse of the mean 
# mu is calculated (i.e., Eq. 2 in manuscript)
BC_h <- 
  as.numeric(sex_specific_mating_system[which(
    sex_specific_mating_system$sex == "female" &
    sex_specific_mating_system$species == "BC"), "mean_annual_no_mates"])
WBC_h <- 
  as.numeric(sex_specific_mating_system[which(
    sex_specific_mating_system$sex == "female" &
    sex_specific_mating_system$species == "WBC"), "mean_annual_no_mates"])

BC_h_sd <- 
  as.numeric(sex_specific_mating_system[which(
    sex_specific_mating_system$sex == "female" &
      sex_specific_mating_system$species == "BC"), "sd_annual_no_mates"])

WBC_h_sd <- 
  as.numeric(sex_specific_mating_system[which(
    sex_specific_mating_system$sex == "female" &
      sex_specific_mating_system$species == "WBC"), "sd_annual_no_mates"])

coucal_mating_system <- 
  data.frame(trait = c("mating_system"),
             species = c("BC", "WBC"),
             mean = c(BC_h, WBC_h),
             sd = c(BC_h_sd, WBC_h_sd))

mating_dat$species <- 
  factor(mating_dat$species ,
         levels = c("BC",
                    "WBC"))

mating_dat_plotting <-
  mating_dat %>% 
  mutate(species_dummy = ifelse(species == "WBC", 1, 0)) %>%
  mutate(species_plot = ifelse(species == "WBC", 2.1, 0.9))

mating_system_plot <-
  ggplot2::ggplot() +
  theme_bw() +
  geom_jitter(aes(y = Nr_partners, x = species_plot, 
                  color = sex),
              data = mating_dat_plotting, 
              width = 0.1, height = 0.1,
              alpha = 0.75, 
              fill = "white", shape = 16) +
  geom_pointrange(data = sex_specific_mating_system, 
    aes(y = mean_annual_no_mates, x = species_plot, 
        ymin = (mean_annual_no_mates-sd_annual_no_mates), 
        ymax = (mean_annual_no_mates+sd_annual_no_mates),
        fill = sex),
    size = 0.8, shape = 21, color = "black") +
  geom_text(aes(y = 0.2, x = species_lab, label = sample_size), 
            data = sex_specific_mating_system,
            size = 3, family = "Verdana") +
  geom_hline(yintercept = 1.5, linetype = "dashed",
             alpha = 0.5, color = "grey20") +
  facet_grid(. ~ sex, labeller = as_labeller(sex_names)) +
  luke_theme +
  theme(legend.position = "none",
        axis.title.x = element_blank(),
        strip.background = element_blank(),
        strip.text = element_text(size = 12, face = "italic")) +
  scale_x_continuous(limits = c(0.5, 2.5), 
                   breaks = c(1, 2), labels = species_names) +
  scale_y_continuous(limits = c(0, 5.5), expand = c(0, 0), 
                     breaks = c(1, 2, 3, 4, 5),
                     labels = c(1, 2, 3, 4, 5)) +
  ylab("Number of unique mates\nper year (± 1 SD)") +
  scale_color_manual(values = sex_pal2) +
  scale_fill_manual(values = sex_pal2) +
  annotate(geom = "text", y = 1.25, x = 1.5,
           label = "monogamy", family = "Verdana", 
           color = "black", size = 3, fontface = 'italic', hjust = 0.5) +
  annotate(geom = "text", y = 1.75, x = 1.5,
           label = "polygamy", family = "Verdana",
           color = "black", size = 3, fontface = 'italic', hjust = 0.5)

Clutch size

Code
# import raw csv data into R
egg_data <-
  read_xls("data/raw/Egg_surv_data_2001_2019_20210524.xls", na = "NA", col_types = "text") %>% 
  dplyr::select(species, nest_ID, year, `Total clutch size`, all_chicks, 
                nst_suc_sho, nest_success, `male_chi+eggs`, `female_chi+eggs`,
                `unkn_sex_eggs+chiks`) %>% 
  filter(species != "CTC") %>% 
  filter(!is.na(nest_success)) %>% 
  dplyr::rename(n_eggs = `Total clutch size`,
                n_chicks = all_chicks,
                male_eggs = `male_chi+eggs`, 
                female_eggs = `female_chi+eggs`,
                unknown_eggs = `unkn_sex_eggs+chiks`) %>% 
  dplyr::mutate(n_eggs = na_if(n_eggs, "?"),
                n_chicks = na_if(n_chicks, "?")) %>% 
  dplyr::mutate(n_eggs = as.numeric(n_eggs),
                n_chicks = as.numeric(n_chicks),
                male_eggs = as.numeric(male_eggs),
                female_eggs = as.numeric(female_eggs),
                unknown_eggs = as.numeric(unknown_eggs)) %>% 
  filter(n_eggs >= n_chicks)

#### clutch size summary ----
coucal_clutch_size <- 
  data.frame(trait = c("clutch_size"),
             species = c("BC", "WBC"),
             mean = c(filter(egg_data, species == "BC") %>% 
                                    summarise(mean_n_eggs = mean(n_eggs, na.rm = TRUE)) %>% 
                                    pull(mean_n_eggs),
                                  filter(egg_data, species == "WBC") %>% 
                                    summarise(mean_n_eggs = mean(n_eggs, na.rm = TRUE)) %>% 
                                    pull(mean_n_eggs)),
             sd = c(filter(egg_data, species == "BC") %>% 
                                  summarise(mean_n_eggs = sd(n_eggs, na.rm = TRUE)) %>% 
                                  pull(mean_n_eggs),
                                filter(egg_data, species == "WBC") %>% 
                                  summarise(mean_n_eggs = sd(n_eggs, na.rm = TRUE)) %>% 
                                  pull(mean_n_eggs)),
             n_nests = c(filter(egg_data, species == "BC") %>% 
                           summarise(n_nests = n_distinct(nest_ID)) %>% 
                           pull(n_nests),
                         filter(egg_data, species == "WBC") %>% 
                           summarise(n_nests = n_distinct(nest_ID)) %>% 
                           pull(n_nests)),
             n_years = c(filter(egg_data, species == "BC") %>% 
                           summarise(n_nests = n_distinct(year)) %>% 
                           pull(n_nests),
                         filter(egg_data, species == "WBC") %>% 
                           summarise(n_nests = n_distinct(year)) %>% 
                           pull(n_nests)))

Egg hatchability

Code
# import raw csv data into R
egg_data <-
  read_xls("data/raw/Egg_surv_data_2001_2019_20210524.xls", 
           na = "NA", col_types = "text") %>% 
  dplyr::select(species, nest_ID, year, `Total clutch size`, all_chicks, 
                nst_suc_sho, nest_success, `male_chi+eggs`, `female_chi+eggs`,
                `unkn_sex_eggs+chiks`) %>% 
  filter(species != "CTC") %>% 
  filter(!is.na(nest_success)) %>% 
  dplyr::rename(n_eggs = `Total clutch size`,
                n_chicks = all_chicks,
                male_eggs = `male_chi+eggs`, 
                female_eggs = `female_chi+eggs`,
                unknown_eggs = `unkn_sex_eggs+chiks`) %>% 
  dplyr::mutate(n_eggs = na_if(n_eggs, "?"),
                n_chicks = na_if(n_chicks, "?")) %>%
  dplyr::mutate(n_eggs = as.numeric(n_eggs),
                n_chicks = as.numeric(n_chicks),
                male_eggs = as.numeric(male_eggs),
                female_eggs = as.numeric(female_eggs),
                unknown_eggs = as.numeric(unknown_eggs)) %>% 
  filter(n_eggs >= n_chicks) %>% 
  arrange(nest_ID)

#### egg hatchability summary ----
egg_surv_data <-
  egg_data %>% 
  dplyr::mutate(n_dead_eggs = n_eggs - n_chicks,
                n_alive_eggs = n_chicks) %>% 
  arrange(nest_ID)

#### sample size summary ----
# egg_surv_data %>% 
#   dplyr::group_by(species) %>% 
#   dplyr::summarise(n_nests = n_distinct(nest_ID))

BC_egg_survival <- 
  lme4::glmer(cbind(n_alive_eggs, n_dead_eggs) ~ 
                1 + 
                (1| nest_ID), 
              data = filter(egg_surv_data, species == "BC"), 
              family = binomial)

WBC_egg_survival <- 
  lme4::glmer(cbind(n_alive_eggs, n_dead_eggs) ~ 
                1 + 
                (1| nest_ID), 
              data = filter(egg_surv_data, species == "WBC"), 
              family = binomial)

# # run tidy bootstrap to obtain model diagnostics
tidy_BC_egg_survival <-
  tidy(BC_egg_survival, conf.int = TRUE, conf.method = "boot", nsim = 1000)

tidy_WBC_egg_survival <-
  tidy(WBC_egg_survival, conf.int = TRUE, conf.method = "boot", nsim = 1000)

coucal_egg_survival <- 
  data.frame(trait = c("egg_survival"),
             species = c("BC", "WBC"),
             mean = c(invlogit(model_parameters(BC_egg_survival)$Coefficient)[1],
                      invlogit(model_parameters(WBC_egg_survival)$Coefficient)[1]),
             CI_low = c(invlogit(model_parameters(BC_egg_survival)$CI_low)[1],
                        invlogit(model_parameters(WBC_egg_survival)$CI_low)[1]),
             CI_high = c(invlogit(model_parameters(BC_egg_survival)$CI_high)[1],
                         invlogit(model_parameters(WBC_egg_survival)$CI_high)[1]),
             n_nests = c(filter(egg_data, species == "BC") %>% 
                           summarise(n_nests = n_distinct(nest_ID)) %>% 
                           pull(n_nests),
                         filter(egg_data, species == "WBC") %>% 
                           summarise(n_nests = n_distinct(nest_ID)) %>% 
                           pull(n_nests)),
             n_years = c(filter(egg_data, species == "BC") %>% 
                           summarise(n_nests = n_distinct(year)) %>% 
                           pull(n_nests),
                         filter(egg_data, species == "WBC") %>% 
                           summarise(n_nests = n_distinct(year)) %>% 
                           pull(n_nests))) %>% 
  mutate(sd = ifelse(!is.na(CI_low), 
                         approx_sd(x1 = CI_low, x2 = CI_high),
                         CI_low))

Hatching sex ratio

Code
egg_data <-
  read_xls("data/raw/Egg_surv_data_2001_2019_20210524.xls", 
           na = "NA", col_types = "text") %>% 
  dplyr::select(species, nest_ID, year, `Total clutch size`, all_chicks, 
                nst_suc_sho, nest_success, `male_chi+eggs`, `female_chi+eggs`,
                `unkn_sex_eggs+chiks`) %>% 
  filter(species != "CTC") %>% 
  filter(!is.na(nest_success)) %>% 
  dplyr::rename(n_eggs = `Total clutch size`,
                n_chicks = all_chicks,
                male_eggs = `male_chi+eggs`, 
                female_eggs = `female_chi+eggs`,
                unknown_eggs = `unkn_sex_eggs+chiks`) %>% 
  dplyr::mutate(n_eggs = na_if(n_eggs, "?"),
                n_chicks = na_if(n_chicks, "?")) %>% 
  dplyr::mutate(n_eggs = as.numeric(n_eggs),
                n_chicks = as.numeric(n_chicks),
                male_eggs = as.numeric(male_eggs),
                female_eggs = as.numeric(female_eggs),
                unknown_eggs = as.numeric(unknown_eggs)) %>% 
  filter(n_eggs >= n_chicks)

#### hatching sex ratio summary ----
egg_HSR_data <- 
  egg_data %>% 
  dplyr::filter(n_eggs == (male_eggs + female_eggs + unknown_eggs)) %>% 
  dplyr::filter(unknown_eggs == 0) %>% 
  mutate(HSR = male_eggs/(male_eggs + female_eggs + unknown_eggs)) %>%
  mutate(species_plot = ifelse(species == "WBC", 2.2, 0.8),
         species_plot2 = ifelse(species == "WBC", 2, 1))

# sample size summary
# egg_HSR_data %>% 
#   dplyr::group_by(species) %>% 
#   dplyr::summarise(n_distinct(nest_ID))

mod_HSR_BC <- 
  lme4::glmer(cbind(male_eggs, female_eggs) ~ 
                1 + 
                (1| nest_ID), 
              data = filter(egg_HSR_data, species == "BC"), 
              family = binomial)

mod_HSR_WBC <- 
  lme4::glmer(cbind(male_eggs, female_eggs) ~ 
                1 + 
                (1| nest_ID), 
              data = filter(egg_HSR_data, species == "WBC"), 
              family = binomial)

# # run tidy bootstrap to obtain model diagnostics
tidy_mod_HSR_BC <-
  tidy(mod_HSR_BC, conf.int = TRUE, conf.method = "boot", nsim = 1000)

tidy_mod_HSR_WBC <-
  tidy(mod_HSR_WBC, conf.int = TRUE, conf.method = "boot", nsim = 1000)

coucal_HSR <- 
  data.frame(trait = "hatching_sex_ratio",
             species = c("BC", "WBC"),
             mean = c(invlogit(model_parameters(mod_HSR_BC)$Coefficient)[1],
                          invlogit(model_parameters(mod_HSR_WBC)$Coefficient)[1]),
             CI_low = c(invlogit(model_parameters(mod_HSR_BC)$CI_low)[1],
                        invlogit(model_parameters(mod_HSR_WBC)$CI_low)[1]),
             CI_high = c(invlogit(model_parameters(mod_HSR_BC)$CI_high)[1],
                         invlogit(model_parameters(mod_HSR_WBC)$CI_high)[1]),
             n_nests = c(filter(egg_HSR_data, species == "BC") %>% 
                           summarise(n_nests = n_distinct(nest_ID)) %>% 
                           pull(n_nests),
                         filter(egg_HSR_data, species == "WBC") %>% 
                           summarise(n_nests = n_distinct(nest_ID)) %>% 
                           pull(n_nests)),
             n_years = c(filter(egg_HSR_data, species == "BC") %>% 
                           summarise(n_nests = n_distinct(year)) %>% 
                           pull(n_nests),
                         filter(egg_HSR_data, species == "WBC") %>% 
                           summarise(n_nests = n_distinct(year)) %>% 
                           pull(n_nests))) %>% 
  mutate(sd = ifelse(!is.na(CI_low), 
                         approx_sd(x1 = CI_low, x2 = CI_high),
                         CI_low),
         species_plot = ifelse(species == "WBC", 1.8, 1.2))

HSR_plot <-
  ggplot2::ggplot() +
  geom_boxplot(data = egg_HSR_data,
               aes(x = species_plot, y = HSR,
                   group = species, fill = species),
               color = "grey50",
               width = 0.05, alpha = 0.5,
               position = position_dodge(width = 0)) +
  geom_errorbar(data = coucal_HSR,
                aes(x = species_plot, ymax = CI_high, ymin = CI_low),
                alpha = 1, color = "black", width = 0.05, lwd = 0.5) +
  geom_point(data = coucal_HSR,
             aes(x = species_plot, y = mean, fill = species),
             lwd = 3, shape = 21, color= "black") +
  geom_jitter(data = egg_HSR_data,
              aes(x = species_plot2, y = HSR,
                  group = species,
                  fill = species, color = species),
              width = 0.02, alpha = 0.2, shape = 19) +
  luke_theme +
  theme(legend.position = "none",
        panel.border = element_blank(),
        strip.background = element_blank(),
        axis.title.x = element_blank(),
        axis.text.x = element_text(12),

        legend.background = element_blank(),
        panel.grid.major.y = element_line(colour = "grey70", size = 0.1),
        plot.title = element_text(hjust = 0.5, face = "italic")) +
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, 0.25)) +
  scale_x_discrete(labels = c("1" = "Black Coucal",
                              "2" = "White-browed Coucal")) +
  ylab(expression(paste("Hatching sex ratio (prop. male)" %+-%  "95% CI", sep = ""))) +
  xlab("Species") +
  scale_color_manual(values = c(brewer.pal(6, "Dark2")[1], brewer.pal(6, "Set1")[1])) +
  scale_fill_manual(values = c(brewer.pal(6, "Dark2")[1], brewer.pal(6, "Set1")[1])) +
  ggtitle("Black Coucal           White-browed\n                                Coucal")

Age at fledging (transition from nestling to groundling)

Code
fledge_dat <- 
  read.csv("data/raw/Coucal_chick_survival_2001-2019_20200129.csv", 
           header = TRUE, stringsAsFactors = FALSE, na.strings = c("", " ", "NA")) %>%
  dplyr::select(species, sex, Fledge_age, nest_ID, year, lab_no, Ring_ID, site) %>% 
  filter(site == "Kapunga" & species != "CTC") %>% 
  mutate(Nst_No = str_replace_all(string = nest_ID, fixed(" "), ""),
         Ring_ID = str_replace_all(string = Ring_ID, fixed(" "), ""),
         lab_no = str_replace_all(string = lab_no, fixed(" "), "")) %>% 
  filter(!is.na(Fledge_age)) %>% 
  mutate(Ind_ID = paste(Nst_No, Ring_ID, lab_no, sep = "_")) %>%
  mutate(sex_plot = ifelse(sex == "M", 2.2, 0.8))
 
#### Modeling (BC) ----
# "Fledge_age" as dependent variable, interaction with sex and species
mod_fledge_age_BC <- 
  lmer(Fledge_age ~ sex + 
         (1 | nest_ID) + (1 | year), 
       data = filter(fledge_dat, species == "BC"))

new_data_BC <- 
  expand.grid(sex = c("F","M"))

predict_BC <- 
  predict(mod_fledge_age_BC, 
          newdata = new_data_BC, 
          re.form = NA, 
          se.fit = TRUE, 
          nsim = 1000)

#### Modeling (WBC) ----
# "Fledge_age" as dependent variable, interaction with sex and species
mod_fledge_age_WBC <- 
  lmer(Fledge_age ~ sex + 
         (1 | nest_ID) + (1 | year), 
       data = filter(fledge_dat, species == "WBC"))

new_data_WBC <- 
  expand.grid(sex = c("F","M"))

predict_WBC <- 
  predict(mod_fledge_age_WBC, 
          newdata = new_data_WBC, 
          re.form = NA, 
          se.fit = TRUE, 
          nsim = 1000)

# tidy up estimates
coucal_fledge_age <- 
  data.frame(trait = c("fledge_age"),
             species = c("BC", "BC", "WBC", "WBC"),
             sex = c("F", "M", "F", "M"),
             mean = c(predict_BC$fit[1],
                      predict_BC$fit[2],
                      predict_WBC$fit[1],
                      predict_WBC$fit[2]),
             CI_low = c(predict_BC$ci.fit[1, 1],
                        predict_BC$ci.fit[1, 2],
                        predict_WBC$ci.fit[1, 1],
                        predict_WBC$ci.fit[1, 2]),
             CI_high = c(predict_BC$ci.fit[2, 1],
                         predict_BC$ci.fit[2, 2],
                         predict_WBC$ci.fit[2, 1],
                         predict_WBC$ci.fit[2, 2]),
             n_inds = c(filter(fledge_dat, species == "BC" & sex == "F") %>% 
                          summarise(n_ = n_distinct(Ind_ID)) %>% 
                          pull(n_),
                        filter(fledge_dat, species == "BC" & sex == "M") %>% 
                          summarise(n_ = n_distinct(Ind_ID)) %>% 
                          pull(n_),
                        filter(fledge_dat, species == "WBC" & sex == "F") %>% 
                          summarise(n_ = n_distinct(Ind_ID)) %>% 
                          pull(n_),
                        filter(fledge_dat, species == "BC" & sex == "M") %>% 
                          summarise(n_ = n_distinct(Ind_ID)) %>% 
                          pull(n_)),
             n_nests = c(filter(fledge_dat, species == "BC" & sex == "F") %>% 
                           summarise(n_ = n_distinct(Nst_No)) %>% 
                           pull(n_),
                         filter(fledge_dat, species == "BC" & sex == "M") %>% 
                           summarise(n_ = n_distinct(Nst_No)) %>% 
                           pull(n_),
                         filter(fledge_dat, species == "WBC" & sex == "F") %>% 
                           summarise(n_ = n_distinct(Nst_No)) %>% 
                           pull(n_),
                         filter(fledge_dat, species == "WBC" & sex == "M") %>% 
                           summarise(n_ = n_distinct(Nst_No)) %>% 
                           pull(n_)),
             n_years = c(filter(fledge_dat, species == "BC" & sex == "F") %>% 
                           summarise(n_ = n_distinct(year)) %>% 
                           pull(n_),
                         filter(fledge_dat, species == "BC" & sex == "M") %>% 
                           summarise(n_ = n_distinct(year)) %>% 
                           pull(n_),
                         filter(fledge_dat, species == "WBC" & sex == "F") %>% 
                           summarise(n_ = n_distinct(year)) %>% 
                           pull(n_),
                         filter(fledge_dat, species == "WBC" & sex == "M") %>% 
                           summarise(n_ = n_distinct(year)) %>% 
                           pull(n_))) %>% 
  mutate(sd = ifelse(!is.na(CI_low), 
                         approx_sd(x1 = CI_low, x2 = CI_high),
                         CI_low))

Black Coucal

Code
#### Modeling (BC) ----
# "Fledge_age" as dependent variable, independent "sex"
mod_fledge_age_BC <- 
  lmer(Fledge_age ~ sex + 
         (1 | nest_ID) + (1 | year), 
       data = filter(fledge_dat, species == "BC"))

# extract model coefficients
mod_fledge_age_BC_coefs <- 
  model_parameters(mod_fledge_age_BC) %>%
  as.data.frame(.)

# run tidy bootstrap to obtain model diagnostics
tidy_mod_fledge_age_BC <-
  tidy(mod_fledge_age_BC, conf.int = TRUE, conf.method = "boot", nsim = 1000)

# run rptR to obtain repeatabilities of random effects
rpt_mod_fledge_age_BC <-
  rpt(Fledge_age ~ sex + 
        (1 | nest_ID) + (1 | year),
      grname = c("nest_ID", "year", "Fixed"),
      data = filter(fledge_dat, species == "BC"),
      datatype = "Gaussian",
      nboot = 1000, npermut = 1000, ratio = TRUE,
      adjusted = TRUE, ncores = 4, parallel = TRUE)

# run partR2 to obtain marginal R2, parameter estimates, and beta weights
R2m_mod_fledge_age_BC <-
  partR2(mod_fledge_age_BC,
         partvars = c("sex"),
         R2_type = "marginal",
         nboot = 1000, CI = 0.95, max_level = 1)

R2c_mod_fledge_age_BC <-
  partR2(mod_fledge_age_BC,
         partvars = c("sex"),
         R2_type = "conditional",
         nboot = 1000, CI = 0.95, max_level = 1)

# save model, tidy, rptR, and partR2 output as a list
stats_mod_fledge_age_BC <-
  list(mod = mod_fledge_age_BC,
       tidy = tidy_mod_fledge_age_BC,
       rptR = rpt_mod_fledge_age_BC,
       partR2m = R2m_mod_fledge_age_BC,
       partR2c = R2c_mod_fledge_age_BC)
Code
#### Table of effect sizes (BC) ----
# Retrieve sample sizes
sample_sizes <-
  fledge_dat %>% 
  filter(species == "BC") %>% 
  summarise(Year = n_distinct(year),
            Individual = n_distinct(Ring_ID),
            Nests = n_distinct(nest_ID))

sample_sizes <- 
  as.data.frame(t(as.data.frame(sample_sizes))) %>%
  rownames_to_column("term") %>% 
  dplyr::rename(estimate = V1) %>% 
  mutate(stat = "n")

# clean model component names
mod_comp_names <- 
  data.frame(comp_name = c("Sex (Male)",
                           "Total Marginal \U1D479\U00B2",
                           "Sex",
                           "Total Conditional \U1D479\U00B2",
                           "Nest",
                           "Year",
                           "Residual",
                           "Nest",
                           "Year",
                           "Residual",
                           "Years",
                           "Individuals",
                           "Observations (i.e., Nests)"))

# Fixed effect sizes (non-standardized)
fixefTable <- 
  stats_mod_fledge_age_BC$tidy %>% 
  dplyr::filter(effect == "fixed") %>% 
  dplyr::select(term, estimate, conf.low, conf.high) %>% 
  as.data.frame() %>% 
  mutate(stat = "fixed")

# Fixed effect sizes (standardized)
fixef_bw_Table <- 
  stats_mod_fledge_age_BC$partR2m$BW %>% 
  # dplyr::select(term, estimate, CI_lower, CI_upper) %>% 
  as.data.frame() %>% 
  mutate(stat = "fixed_bw") %>% 
  dplyr::rename(conf.low = CI_lower,
                conf.high = CI_upper)

# Semi-partial R2 estimates
fixef_bw_Table <-
  model_parameters(stats_mod_fledge_age_BC$mod, standardize = "refit") %>%
  dplyr::select(Parameter, Coefficient, CI_low, CI_high) %>% 
  as.data.frame() %>% 
  mutate(stat = "fixed_bw") %>% 
  dplyr::rename(conf.low = CI_low,
                conf.high = CI_high,
                term = Parameter,
                estimate = Coefficient) %>% 
  filter(term == "sexM")

R2Table <- 
  bind_rows(stats_mod_fledge_age_BC$partR2m$R2,
            stats_mod_fledge_age_BC$partR2c$R2[1,]) %>% 
  dplyr::select(term, estimate, CI_lower, CI_upper) %>% 
  as.data.frame() %>% 
  mutate(stat = "partR2") %>% 
  dplyr::rename(conf.low = CI_lower,
                conf.high = CI_upper)

# Random effects variances
ranefTable <- 
  stats_mod_fledge_age_BC$tidy %>% 
  dplyr::filter(effect == "ran_pars") %>% 
  dplyr::select(group, estimate, conf.low, conf.high) %>% 
  as.data.frame() %>% 
  mutate(stat = "rand") %>% 
  dplyr::rename(term = group) %>% 
  mutate(estimate = estimate^2,
         conf.high = conf.high^2,
         conf.low = conf.low^2)

# Adjusted repeatabilities
coefRptTable <- 
  stats_mod_fledge_age_BC$rptR$R_boot %>% 
  dplyr::select(-Fixed) %>% 
  mutate(residual = 1 - rowSums(.)) %>% 
  apply(., 2, 
        function(x) c(mean (x), quantile (x, prob = c(0.025, 0.975)))) %>% 
  t() %>% 
  as.data.frame() %>% 
  rownames_to_column("term") %>% 
  dplyr::rename(estimate = V1,
                conf.low = `2.5%`,
                conf.high = `97.5%`) %>% 
  mutate(stat = "RptR")

# Store all parameters into a single table and clean it up
allCoefs_mod <- 
  bind_rows(fixef_bw_Table,
            R2Table,
            ranefTable, 
            coefRptTable, 
            sample_sizes) %>% 
  bind_cols(.,
            mod_comp_names) %>%
  mutate(coefString = ifelse(!is.na(conf.low),
                             paste0("[", 
                                    round(conf.low, 2), ", ", 
                                    round(conf.high, 2), "]"),
                             NA),
         effect = c(rep("Fixed effects \U1D6FD (standardized)", nrow(fixef_bw_Table)),
                    rep("Partitioned \U1D479\U00B2", nrow(R2Table)),
                    rep("Random effects \U1D70E\U00B2", nrow(ranefTable)),
                    rep("Adjusted repeatability \U1D45F", nrow(coefRptTable)),
                    rep("Sample sizes \U1D45B", nrow(sample_sizes)))) %>%
  dplyr::select(effect, everything())

# draw gt table
table_mod_fledge_age_BC <- 
  allCoefs_mod %>% 
  dplyr::select(effect, comp_name, estimate, coefString) %>% 
  gt(rowname_col = "row",
     groupname_col = "effect") %>% 
  cols_width(vars(comp_name) ~ pct(50),
             vars(estimate) ~ pct(20),
             vars(coefString) ~ pct(30)) %>% 
  cols_label(comp_name = html("<i>Black Coucals <br> age at leaving nest</i>"),
             estimate = "Mean estimate",
             coefString = "95% confidence interval") %>% 
  fmt_number(columns = vars(estimate),
             rows = 1:10,
             decimals = 2,
             use_seps = FALSE) %>% 
  fmt_number(columns = vars(estimate),
             rows = 11:13,
             decimals = 0,
             use_seps = FALSE) %>% 
  fmt_missing(columns = 1:4,
              missing_text = "") %>% 
  cols_align(align = "left",
             columns = vars(comp_name)) %>% 
  tab_options(row_group.font.weight = "bold",
              row_group.background.color = brewer.pal(9,"Greys")[3],
              table.font.size = 12,
              data_row.padding = 3,
              row_group.padding = 4,
              summary_row.padding = 2,
              column_labels.font.size = 14,
              row_group.font.size = 12,
              table.width = pct(80))

table_mod_fledge_age_BC
Black Coucals
age at leaving nest
Mean estimate 95% confidence interval
Fixed effects 𝛽 (standardized)
Sex (Male) −0.31 [-0.47, -0.15]
Partitioned 𝑹²
Total Marginal 𝑹² 0.02 [0.01, 0.05]
Sex 0.02 [0.01, 0.05]
Total Conditional 𝑹² 0.46 [0.35, 0.56]
Random effects 𝜎²
Nest 1.06 [0.7, 1.47]
Year 0.06 [0, 0.29]
Residual 1.38 [1.15, 1.6]
Adjusted repeatability 𝑟
Nest 0.42 [0.32, 0.52]
Year 0.03 [0, 0.11]
Residual 0.55 [0.45, 0.66]
Sample sizes 𝑛
Years 13
Individuals 216
Observations (i.e., Nests) 133

White-browed coucal

Code
#### Modeling (WBC) ----
# "Fledge_age" as dependent variable, interaction with sex and species
mod_fledge_age_WBC <- 
  lmer(Fledge_age ~ sex + 
         (1 | nest_ID) + (1 | year), 
       data = filter(fledge_dat, species == "WBC"))

# extract model coefficients
mod_fledge_age_WBC_coefs <- 
  model_parameters(mod_fledge_age_WBC) %>%
  as.data.frame(.)

# run tidy bootstrap to obtain model diagnostics
tidy_mod_fledge_age_WBC <-
  tidy(mod_fledge_age_WBC, conf.int = TRUE, conf.method = "boot", nsim = 1000)

# run rptR to obtain repeatabilities of random effects
rpt_mod_fledge_age_WBC <-
  rpt(Fledge_age ~ sex + 
        (1 | nest_ID) + (1 | year),
      grname = c("nest_ID", "year", "Fixed"),
      data = filter(fledge_dat, species == "WBC"),
      datatype = "Gaussian",
      nboot = 1000, npermut = 1000, ratio = TRUE,
      adjusted = TRUE, ncores = 4, parallel = TRUE)

# run partR2 to obtain marginal R2, parameter estimates, and beta weights
R2m_mod_fledge_age_WBC <-
  partR2(mod_fledge_age_WBC,
         partvars = c("sex"),
         R2_type = "marginal",
         nboot = 1000, CI = 0.95, max_level = 1)

R2c_mod_fledge_age_WBC <-
  partR2(mod_fledge_age_WBC,
         partvars = c("sex"),
         R2_type = "conditional",
         nboot = 1000, CI = 0.95, max_level = 1)

# save model, tidy, rptR, and partR2 output as a list
stats_mod_fledge_age_WBC <-
  list(mod = mod_fledge_age_WBC,
       tidy = tidy_mod_fledge_age_WBC,
       rptR = rpt_mod_fledge_age_WBC,
       partR2m = R2m_mod_fledge_age_WBC,
       partR2c = R2c_mod_fledge_age_WBC)
Code
#### Table of effect sizes (WBC) ----
# Retrieve sample sizes
sample_sizes <-
  fledge_dat %>% 
  filter(species == "WBC") %>% 
  summarise(Year = n_distinct(year),
            Individual = n_distinct(Ring_ID),
            Nests = n_distinct(nest_ID))

sample_sizes <- 
  as.data.frame(t(as.data.frame(sample_sizes))) %>%
  rownames_to_column("term") %>% 
  dplyr::rename(estimate = V1) %>% 
  mutate(stat = "n")

# clean model component names
mod_comp_names <- 
  data.frame(comp_name = c("Sex (Male)",
                           "Total Marginal \U1D479\U00B2",
                           "Sex",
                           "Total Conditional \U1D479\U00B2",
                           "Nest",
                           "Year",
                           "Residual",
                           "Nest",
                           "Year",
                           "Residual",
                           "Years",
                           "Individuals",
                           "Observations (i.e., Nests)"))

# Fixed effect sizes (non-standardized)
fixefTable <- 
  stats_mod_fledge_age_WBC$tidy %>% 
  dplyr::filter(effect == "fixed") %>% 
  dplyr::select(term, estimate, conf.low, conf.high) %>% 
  as.data.frame() %>% 
  mutate(stat = "fixed")

# Fixed effect sizes (standardized)
fixef_bw_Table <- 
  stats_mod_fledge_age_WBC$partR2m$BW %>% 
  # dplyr::select(term, estimate, CI_lower, CI_upper) %>% 
  as.data.frame() %>% 
  mutate(stat = "fixed_bw") %>% 
  dplyr::rename(conf.low = CI_lower,
                conf.high = CI_upper)

# Semi-partial R2 estimates
fixef_bw_Table <-
  model_parameters(stats_mod_fledge_age_WBC$mod, standardize = "refit") %>%
  dplyr::select(Parameter, Coefficient, CI_low, CI_high) %>% 
  as.data.frame() %>% 
  mutate(stat = "fixed_bw") %>% 
  dplyr::rename(conf.low = CI_low,
                conf.high = CI_high,
                term = Parameter,
                estimate = Coefficient) %>% 
  filter(term == "sexM")

R2Table <- 
  bind_rows(stats_mod_fledge_age_WBC$partR2m$R2,
            stats_mod_fledge_age_WBC$partR2c$R2[1,]) %>% 
  dplyr::select(term, estimate, CI_lower, CI_upper) %>% 
  as.data.frame() %>% 
  mutate(stat = "partR2") %>% 
  dplyr::rename(conf.low = CI_lower,
                conf.high = CI_upper)

# Random effects variances
ranefTable <- 
  stats_mod_fledge_age_WBC$tidy %>% 
  dplyr::filter(effect == "ran_pars") %>% 
  dplyr::select(group, estimate, conf.low, conf.high) %>% 
  as.data.frame() %>% 
  mutate(stat = "rand") %>% 
  dplyr::rename(term = group) %>% 
  mutate(estimate = estimate^2,
         conf.high = conf.high^2,
         conf.low = conf.low^2)

# Adjusted repeatabilities
coefRptTable <- 
  stats_mod_fledge_age_WBC$rptR$R_boot %>% 
  dplyr::select(-Fixed) %>% 
  mutate(residual = 1 - rowSums(.)) %>% 
  apply(., 2, 
        function(x) c(mean (x), quantile (x, prob = c(0.025, 0.975)))) %>% 
  t() %>% 
  as.data.frame() %>% 
  rownames_to_column("term") %>% 
  dplyr::rename(estimate = V1,
                conf.low = `2.5%`,
                conf.high = `97.5%`) %>% 
  mutate(stat = "RptR")

# Store all parameters into a single table and clean it up
allCoefs_mod <- 
  bind_rows(fixef_bw_Table,
            R2Table,
            ranefTable, 
            coefRptTable, 
            sample_sizes) %>% 
  bind_cols(.,
            mod_comp_names) %>%
  mutate(coefString = ifelse(!is.na(conf.low),
                             paste0("[", 
                                    round(conf.low, 2), ", ", 
                                    round(conf.high, 2), "]"),
                             NA),
         effect = c(rep("Fixed effects \U1D6FD (standardized)", nrow(fixef_bw_Table)),
                    rep("Partitioned \U1D479\U00B2", nrow(R2Table)),
                    rep("Random effects \U1D70E\U00B2", nrow(ranefTable)),
                    rep("Adjusted repeatability \U1D45F", nrow(coefRptTable)),
                    rep("Sample sizes \U1D45B", nrow(sample_sizes)))) %>%
  dplyr::select(effect, everything())

# draw gt table
table_mod_fledge_age_WBC <- 
  allCoefs_mod %>% 
  dplyr::select(effect, comp_name, estimate, coefString) %>% 
  gt(rowname_col = "row",
     groupname_col = "effect") %>% 
  cols_width(vars(comp_name) ~ pct(50),
             vars(estimate) ~ pct(20),
             vars(coefString) ~ pct(30)) %>%
  cols_label(comp_name = html("<i>White-browed Coucals <br> age at leaving nest</i>"),
             estimate = "Mean estimate",
             coefString = "95% confidence interval") %>% 
  fmt_number(columns = vars(estimate),
             rows = 1:10,
             decimals = 2,
             use_seps = FALSE) %>% 
  fmt_number(columns = vars(estimate),
             rows = 11:13,
             decimals = 0,
             use_seps = FALSE) %>% 
  fmt_missing(columns = 1:4,
              missing_text = "") %>% 
  cols_align(align = "left",
             columns = vars(comp_name)) %>% 
  tab_options(row_group.font.weight = "bold",
              row_group.background.color = brewer.pal(9,"Greys")[3],
              table.font.size = 12,
              data_row.padding = 3,
              row_group.padding = 4,
              summary_row.padding = 2,
              column_labels.font.size = 14,
              row_group.font.size = 12,
              table.width = pct(80))
White-browed Coucals
age at leaving nest
Mean estimate 95% confidence interval
Fixed effects 𝛽 (standardized)
Sex (Male) −0.01 [-0.19, 0.16]
Partitioned 𝑹²
Total Marginal 𝑹² 0.00 [0, 0.01]
Sex 0.00 [0, 0.01]
Total Conditional 𝑹² 0.58 [0.46, 0.67]
Random effects 𝜎²
Nest 2.13 [1.41, 2.87]
Year 0.00 [0, 0.33]
Residual 1.56 [1.29, 1.85]
Adjusted repeatability 𝑟
Nest 0.56 [0.45, 0.66]
Year 0.01 [0, 0.08]
Residual 0.43 [0.33, 0.53]
Sample sizes 𝑛
Years 12
Individuals 279
Observations (i.e., Nests) 109
Code
#### Forest plot of model predictions ----
# extract model coefficients
mod_fledge_age_BC_coefs <- 
  model_parameters(mod_fledge_age_BC) %>%
  as.data.frame(.)

mod_fledge_age_WBC_coefs <- 
  model_parameters(mod_fledge_age_BC) %>%
  as.data.frame(.)

mod_fledge_age_BC_fits <- 
  as.data.frame(effect(term = "sex", mod = stats_mod_fledge_age_BC$mod, 
                       xlevels = list(sex = c("M", "F")))) %>%
  mutate(sex_plot = ifelse(sex == "M", 1.8, 1.2),
         species = "BC")

mod_fledge_age_WBC_fits <- 
  as.data.frame(effect(term = "sex", mod = stats_mod_fledge_age_WBC$mod, 
                       xlevels = list(sex = c("M", "F")))) %>%
  mutate(sex_plot = ifelse(sex == "M", 1.8, 1.2),
         species = "WBC")

mod_fledge_age_fits <- 
  bind_rows(mod_fledge_age_BC_fits, 
            mod_fledge_age_WBC_fits) %>% 
  `rownames<-`( NULL )

# plot_palette_recruit <- brewer.pal(6, "Dark2")[c(2,3)]

# mod_fledge_age_fits[1, 2] - mod_fledge_age_fits[2, 2] 
# mod_fledge_age_fits[1, 4] - mod_fledge_age_fits[2, 4]
# mod_fledge_age_fits[1, 5] - mod_fledge_age_fits[2, 5] 

fledge_age_plot <-
  ggplot2::ggplot() + 
  geom_boxplot(data = fledge_dat,
               aes(x = sex_plot, y = Fledge_age,
                   group = sex, fill = sex),
               color = "grey50",
               width = 0.05, alpha = 0.5,
               position = position_dodge(width = 0)) +
  geom_errorbar(data = mod_fledge_age_fits, 
                aes(x = sex_plot, ymax = upper, ymin = lower),
                alpha = 1, color = "black", width = 0.05, lwd = 0.5) + 
  geom_point(data = mod_fledge_age_fits, 
             aes(x = sex_plot, y = fit, fill = sex),
             lwd = 3, shape = 21, color= "black") +
  geom_jitter(data = fledge_dat, 
              aes(x = sex, y = Fledge_age, 
                  group = sex, 
                  fill = sex, color = sex), 
              width = 0.02, alpha = 0.2, shape = 19) +
  coord_flip() +
  luke_theme +
  theme(legend.position = "none",
        panel.border = element_blank(),
        strip.background = element_blank(),
        strip.text = element_text(size = 12, face = "italic"),
        axis.title.x = element_text(size = 12),
        # axis.title.y = element_blank(),
        # axis.text.y = element_blank(),
        # axis.ticks.x = element_blank(),
        axis.ticks.y = element_blank(),
        legend.background = element_blank(),
        # panel.grid = element_blank(),
        panel.grid.major.x = element_line(colour = "grey70", size = 0.1),
        axis.title.y = element_blank()) +
  scale_y_continuous(limits = c(7.5, 20.5), breaks = seq(8, 20, 2)) +#c(-60, -30, 0, 30, 60)) +
  scale_x_discrete(labels = c("M" = "Male",
                              "F" = "Female")) +
  ylab(expression(paste("Age at leaving nest (days)" %+-%  "95% CI", sep = ""))) +
  # xlab("Origin") +
  scale_color_manual(values = sex_pal2) +
  scale_fill_manual(values = sex_pal2) +
  facet_grid(species ~ ., labeller = as_labeller(species_names))# +
# annotate(geom = "text", x = 0.5, y = 58,
#          label = "First nests of the season",
#          color = "black", size = 3, fontface = 'italic', hjust = 0)

fledge_age_plot

Age at first flight (transition from groundling to fledgling)

Here we assess sex-specific age at first flight while controlling for non-independence among siblings, age at fledging, and mass at fledging.

Code
flight_dat <- 
  read.csv("data/raw/coucals_first_flights_data_per_individual2020.csv", 
           header = TRUE, stringsAsFactors = FALSE, na.strings = c("", " ", "NA")) %>% 
  mutate(Nst_No = str_replace_all(string = Nst_No, fixed(" "), ""),
         Ind_ID = str_replace_all(string = Ind_ID, fixed(" "), "")) %>%
  dplyr::rename(sex = Sex,
                nest_ID = Nst_No,
                year = Year,
                species = Spp,
                flight_age = days_since_fledging,
                Ring_ID = Ind_ID) %>% 
  mutate(sex_plot = ifelse(sex == "Male", 2.2, 0.8))

Black Coucals

Code
#### Modeling (BC) ----
# "flight_age" as dependent variable, 
# "sex", "fledge_age", and "fledge_mass" as independent
mod_flight_age_BC <- 
  lmer(flight_age ~ sex + Fledge_age + Fledge_mass +
         (1 | nest_ID) + (1 | year), 
       data = filter(flight_dat, species == "BC"))

# extract model coefficients
mod_flight_age_BC_coefs <- 
  model_parameters(mod_flight_age_BC) %>%
  as.data.frame(.)

# run tidy bootstrap to obtain model diagnostics
tidy_mod_flight_age_BC <-
  tidy(mod_flight_age_BC, conf.int = TRUE, conf.method = "boot", nsim = 1000)

# run rptR to obtain repeatabilities of random effects
rpt_mod_flight_age_BC <-
  rpt(flight_age ~ sex + Fledge_age + Fledge_mass + 
        (1 | nest_ID) + (1 | year),
      grname = c("nest_ID", "year", "Fixed"),
      data = filter(flight_dat, species == "BC"),
      datatype = "Gaussian",
      nboot = 1000, npermut = 1000, ratio = TRUE,
      adjusted = TRUE, ncores = 4, parallel = TRUE)

# run partR2 to obtain marginal R2, parameter estimates, and beta weights
R2m_mod_flight_age_BC <-
  partR2(mod_flight_age_BC,
         partvars = c("sex", "Fledge_age", "Fledge_mass"),
         R2_type = "marginal",
         nboot = 1000, CI = 0.95, max_level = 1)

R2c_mod_flight_age_BC <-
  partR2(mod_flight_age_BC,
         partvars = c("sex", "Fledge_age", "Fledge_mass"),
         R2_type = "conditional",
         nboot = 1000, CI = 0.95, max_level = 1)

# save model, tidy, rptR, and partR2 output as a list
stats_mod_flight_age_BC <-
  list(mod = mod_flight_age_BC,
       tidy = tidy_mod_flight_age_BC,
       rptR = rpt_mod_flight_age_BC,
       partR2m = R2m_mod_flight_age_BC,
       partR2c = R2c_mod_flight_age_BC)
Code
#### Table of effect sizes (BC) ----
# Retrieve sample sizes
sample_sizes <-
  flight_dat %>% 
  filter(species == "BC") %>% 
  summarise(Year = n_distinct(year),
            Individual = n_distinct(Ring_ID),
            Nests = n_distinct(nest_ID))

sample_sizes <- 
  as.data.frame(t(as.data.frame(sample_sizes))) %>%
  rownames_to_column("term") %>% 
  dplyr::rename(estimate = V1) %>% 
  mutate(stat = "n")

# clean model component names
mod_comp_names <- 
  data.frame(comp_name = c("Sex (Male)",
                           "Age at nest departure",
                           "Mass at nest departure",
                           "Total Marginal \U1D479\U00B2",
                           "Sex",
                           "Age at nest departure",
                           "Mass at nest departure",
                           "Total Conditional \U1D479\U00B2",
                           "Nest",
                           "Year",
                           "Residual",
                           "Nest",
                           "Year",
                           "Residual",
                           "Years",
                           "Individuals",
                           "Observations (i.e., Nests)"))

# Fixed effect sizes (non-standardized)
fixefTable <- 
  stats_mod_flight_age_BC$tidy %>% 
  dplyr::filter(effect == "fixed") %>% 
  dplyr::select(term, estimate, conf.low, conf.high) %>% 
  as.data.frame() %>% 
  mutate(stat = "fixed")

# Fixed effect sizes (standardized)
fixef_bw_Table <- 
  stats_mod_flight_age_BC$partR2m$BW %>% 
  as.data.frame() %>% 
  mutate(stat = "fixed_bw") %>% 
  dplyr::rename(conf.low = CI_lower,
                conf.high = CI_upper)

# Semi-partial R2 estimates
fixef_bw_Table <-
  model_parameters(stats_mod_flight_age_BC$mod, standardize = "refit") %>%
  dplyr::select(Parameter, Coefficient, CI_low, CI_high) %>% 
  as.data.frame() %>% 
  mutate(stat = "fixed_bw") %>% 
  dplyr::rename(conf.low = CI_low,
                conf.high = CI_high,
                term = Parameter,
                estimate = Coefficient) %>% 
  filter(term %in% c("sexMale", "Fledge_age", "Fledge_mass"))

R2Table <- 
  bind_rows(stats_mod_flight_age_BC$partR2m$R2,
            stats_mod_flight_age_BC$partR2c$R2[1,]) %>% 
  dplyr::select(term, estimate, CI_lower, CI_upper) %>% 
  as.data.frame() %>% 
  mutate(stat = "partR2") %>% 
  dplyr::rename(conf.low = CI_lower,
                conf.high = CI_upper)

# Random effects variances
ranefTable <- 
  stats_mod_flight_age_BC$tidy %>% 
  dplyr::filter(effect == "ran_pars") %>% 
  dplyr::select(group, estimate, conf.low, conf.high) %>% 
  as.data.frame() %>% 
  mutate(stat = "rand") %>% 
  dplyr::rename(term = group) %>% 
  mutate(estimate = estimate^2,
         conf.high = conf.high^2,
         conf.low = conf.low^2)

# Adjusted repeatabilities
coefRptTable <- 
  stats_mod_flight_age_BC$rptR$R_boot %>% 
  dplyr::select(-Fixed) %>% 
  mutate(residual = 1 - rowSums(.)) %>% 
  apply(., 2, 
        function(x) c(mean (x), quantile (x, prob = c(0.025, 0.975)))) %>% 
  t() %>% 
  as.data.frame() %>% 
  rownames_to_column("term") %>% 
  dplyr::rename(estimate = V1,
                conf.low = `2.5%`,
                conf.high = `97.5%`) %>% 
  mutate(stat = "RptR")

# Store all parameters into a single table and clean it up
allCoefs_mod <- 
  bind_rows(fixef_bw_Table,
            R2Table,
            ranefTable, 
            coefRptTable, 
            sample_sizes) %>% 
  bind_cols(.,
            mod_comp_names) %>%
  mutate(coefString = ifelse(!is.na(conf.low),
                             paste0("[", 
                                    round(conf.low, 2), ", ", 
                                    round(conf.high, 2), "]"),
                             NA),
         effect = c(rep("Fixed effects \U1D6FD (standardized)", nrow(fixef_bw_Table)),
                    rep("Partitioned \U1D479\U00B2", nrow(R2Table)),
                    rep("Random effects \U1D70E\U00B2", nrow(ranefTable)),
                    rep("Adjusted repeatability \U1D45F", nrow(coefRptTable)),
                    rep("Sample sizes \U1D45B", nrow(sample_sizes)))) %>%
  dplyr::select(effect, everything())

# draw gt table
table_mod_flight_age_BC <- 
  allCoefs_mod %>% 
  dplyr::select(effect, comp_name, estimate, coefString) %>% 
  gt(rowname_col = "row",
     groupname_col = "effect") %>% 
  cols_width(vars(comp_name) ~ pct(50),
             vars(estimate) ~ pct(20),
             vars(coefString) ~ pct(30)) %>% 
  cols_label(comp_name = html("<i>Black Coucals <br> groundling duration</i>"),
             estimate = "Mean estimate",
             coefString = "95% confidence interval") %>% 
  fmt_number(columns = vars(estimate),
             rows = 1:14,
             decimals = 2,
             use_seps = FALSE) %>% 
  fmt_number(columns = vars(estimate),
             rows = 15:17,
             decimals = 0,
             use_seps = FALSE) %>% 
  fmt_missing(columns = 1:4,
              missing_text = "") %>% 
  cols_align(align = "left",
             columns = vars(comp_name)) %>% 
  tab_options(row_group.font.weight = "bold",
              row_group.background.color = brewer.pal(9,"Greys")[3],
              table.font.size = 12,
              data_row.padding = 3,
              row_group.padding = 4,
              summary_row.padding = 2,
              column_labels.font.size = 14,
              row_group.font.size = 12,
              table.width = pct(80))

table_mod_flight_age_BC
Black Coucals
groundling duration
Mean estimate 95% confidence interval
Fixed effects 𝛽 (standardized)
Sex (Male) −0.01 [-0.42, 0.41]
Age at nest departure 0.04 [-0.17, 0.26]
Mass at nest departure 0.14 [-0.1, 0.38]
Partitioned 𝑹²
Total Marginal 𝑹² 0.02 [0, 0.15]
Sex 0.00 [0, 0.12]
Age at nest departure 0.00 [0, 0.12]
Mass at nest departure 0.01 [0, 0.13]
Total Conditional 𝑹² 0.54 [0.29, 0.75]
Random effects 𝜎²
Nest 20.74 [5.6, 37.38]
Year 1.57 [0, 11.35]
Residual 19.90 [11.73, 30.89]
Adjusted repeatability 𝑟
Nest 0.46 [0.18, 0.7]
Year 0.05 [0, 0.23]
Residual 0.49 [0.27, 0.78]
Sample sizes 𝑛
Years 5
Individuals 85
Observations (i.e., Nests) 45

White-browed Coucals

Code
# "flight_age" as dependent variable, 
# "sex", "fledge_age", and "fledge_mass" as independent
mod_flight_age_WBC <-
  lmer(flight_age ~ sex + Fledge_age + Fledge_mass +
         (1 | nest_ID) + (1 | year), 
       data = filter(flight_dat, species == "WBC"))

# extract model coefficients
mod_flight_age_WBC_coefs <- 
  model_parameters(mod_flight_age_WBC) %>%
  as.data.frame(.)

# run tidy bootstrap to obtain model diagnostics
tidy_mod_flight_age_WBC <-
  tidy(mod_flight_age_WBC, conf.int = TRUE, conf.method = "boot", nsim = 1000)

# run rptR to obtain repeatabilities of random effects
rpt_mod_flight_age_WBC <-
  rpt(flight_age ~ sex + Fledge_age + Fledge_mass + 
        (1 | nest_ID) + (1 | year),
      grname = c("nest_ID", "year", "Fixed"),
      data = filter(flight_dat, species == "WBC"),
      datatype = "Gaussian",
      nboot = 1000, npermut = 1000, ratio = TRUE,
      adjusted = TRUE, ncores = 4, parallel = TRUE)

# run partR2 to obtain marginal R2, parameter estimates, and beta weights
R2m_mod_flight_age_WBC <-
  partR2(mod_flight_age_WBC,
         partvars = c("sex", "Fledge_age", "Fledge_mass"),
         R2_type = "marginal",
         nboot = 1000, CI = 0.95, max_level = 1)

R2c_mod_flight_age_WBC <-
  partR2(mod_flight_age_WBC,
         partvars = c("sex", "Fledge_age", "Fledge_mass"),
         R2_type = "conditional",
         nboot = 1000, CI = 0.95, max_level = 1)

# save model, tidy, rptR, and partR2 output as a list
stats_mod_flight_age_WBC <-
  list(mod = mod_flight_age_WBC,
       tidy = tidy_mod_flight_age_WBC,
       rptR = rpt_mod_flight_age_WBC,
       partR2m = R2m_mod_flight_age_WBC,
       partR2c = R2c_mod_flight_age_WBC)
Code
#### Table of effect sizes (BC) ----
# Retrieve sample sizes
sample_sizes <-
  flight_dat %>% 
  filter(species == "WBC") %>% 
  summarise(Year = n_distinct(year),
            Individual = n_distinct(Ring_ID),
            Nests = n_distinct(nest_ID))

sample_sizes <- 
  as.data.frame(t(as.data.frame(sample_sizes))) %>%
  rownames_to_column("term") %>% 
  dplyr::rename(estimate = V1) %>% 
  mutate(stat = "n")

# clean model component names
mod_comp_names <- 
  data.frame(comp_name = c("Sex (Male)",
                           "Age at nest departure",
                           "Mass at nest departure",
                           "Total Marginal \U1D479\U00B2",
                           "Sex",
                           "Age at nest departure",
                           "Mass at nest departure",
                           "Total Conditional \U1D479\U00B2",
                           "Nest",
                           "Year",
                           "Residual",
                           "Nest",
                           "Year",
                           "Residual",
                           "Years",
                           "Individuals",
                           "Observations (i.e., Nests)"))

# Fixed effect sizes (non-standardized)
fixefTable <- 
  stats_mod_flight_age_WBC$tidy %>% 
  dplyr::filter(effect == "fixed") %>% 
  dplyr::select(term, estimate, conf.low, conf.high) %>% 
  as.data.frame() %>% 
  mutate(stat = "fixed")

# Fixed effect sizes (standardized)
fixef_bw_Table <- 
  stats_mod_flight_age_WBC$partR2m$BW %>% 
  # dplyr::select(term, estimate, CI_lower, CI_upper) %>% 
  as.data.frame() %>% 
  mutate(stat = "fixed_bw") %>% 
  dplyr::rename(conf.low = CI_lower,
                conf.high = CI_upper)

# Semi-partial R2 estimates
fixef_bw_Table <-
  model_parameters(stats_mod_flight_age_WBC$mod, standardize = "refit") %>%
  dplyr::select(Parameter, Coefficient, CI_low, CI_high) %>% 
  as.data.frame() %>% 
  mutate(stat = "fixed_bw") %>% 
  dplyr::rename(conf.low = CI_low,
                conf.high = CI_high,
                term = Parameter,
                estimate = Coefficient) %>% 
  filter(term %in% c("sexMale", "Fledge_age", "Fledge_mass"))

R2Table <- 
  bind_rows(stats_mod_flight_age_WBC$partR2m$R2,
            stats_mod_flight_age_WBC$partR2c$R2[1,]) %>% 
  dplyr::select(term, estimate, CI_lower, CI_upper) %>% 
  as.data.frame() %>% 
  mutate(stat = "partR2") %>% 
  dplyr::rename(conf.low = CI_lower,
                conf.high = CI_upper)

# Random effects variances
ranefTable <- 
  stats_mod_flight_age_WBC$tidy %>% 
  dplyr::filter(effect == "ran_pars") %>% 
  dplyr::select(group, estimate, conf.low, conf.high) %>% 
  as.data.frame() %>% 
  mutate(stat = "rand") %>% 
  dplyr::rename(term = group) %>% 
  mutate(estimate = estimate^2,
         conf.high = conf.high^2,
         conf.low = conf.low^2)

# Adjusted repeatabilities
coefRptTable <- 
  stats_mod_flight_age_WBC$rptR$R_boot %>% 
  dplyr::select(-Fixed) %>% 
  mutate(residual = 1 - rowSums(.)) %>% 
  apply(., 2, 
        function(x) c(mean (x), quantile (x, prob = c(0.025, 0.975)))) %>% 
  t() %>% 
  as.data.frame() %>% 
  rownames_to_column("term") %>% 
  dplyr::rename(estimate = V1,
                conf.low = `2.5%`,
                conf.high = `97.5%`) %>% 
  mutate(stat = "RptR")

# Store all parameters into a single table and clean it up
allCoefs_mod <- 
  bind_rows(fixef_bw_Table,
            R2Table,
            ranefTable, 
            coefRptTable, 
            sample_sizes) %>% 
  bind_cols(.,
            mod_comp_names) %>%
  mutate(coefString = ifelse(!is.na(conf.low),
                             paste0("[", 
                                    round(conf.low, 2), ", ", 
                                    round(conf.high, 2), "]"),
                             NA),
         effect = c(rep("Fixed effects \U1D6FD (standardized)", nrow(fixef_bw_Table)),
                    rep("Partitioned \U1D479\U00B2", nrow(R2Table)),
                    rep("Random effects \U1D70E\U00B2", nrow(ranefTable)),
                    rep("Adjusted repeatability \U1D45F", nrow(coefRptTable)),
                    rep("Sample sizes \U1D45B", nrow(sample_sizes)))) %>%
  dplyr::select(effect, everything())

# draw gt table
table_mod_flight_age_WBC <- 
  allCoefs_mod %>% 
  dplyr::select(effect, comp_name, estimate, coefString) %>% 
  gt(rowname_col = "row",
     groupname_col = "effect") %>% 
  cols_width(vars(comp_name) ~ pct(50),
             vars(estimate) ~ pct(20),
             vars(coefString) ~ pct(30)) %>% 
  cols_label(comp_name = html("<i>White-browed Coucals <br> groundling duration</i>"),
             estimate = "Mean estimate",
             coefString = "95% confidence interval") %>% 
  fmt_number(columns = vars(estimate),
             rows = 1:14,
             decimals = 2,
             use_seps = FALSE) %>% 
  fmt_number(columns = vars(estimate),
             rows = 15:17,
             decimals = 0,
             use_seps = FALSE) %>% 
  fmt_missing(columns = 1:4,
              missing_text = "") %>% 
  cols_align(align = "left",
             columns = vars(comp_name)) %>% 
  tab_options(row_group.font.weight = "bold",
              row_group.background.color = brewer.pal(9,"Greys")[3],
              table.font.size = 12,
              data_row.padding = 3,
              row_group.padding = 4,
              summary_row.padding = 2,
              column_labels.font.size = 14,
              row_group.font.size = 12,
              table.width = pct(80))

table_mod_flight_age_WBC
White-browed Coucals
groundling duration
Mean estimate 95% confidence interval
Fixed effects 𝛽 (standardized)
Sex (Male) −0.04 [-0.5, 0.42]
Age at nest departure −0.04 [-0.32, 0.25]
Mass at nest departure 0.17 [-0.1, 0.45]
Partitioned 𝑹²
Total Marginal 𝑹² 0.03 [0.01, 0.22]
Sex 0.00 [0, 0.18]
Age at nest departure 0.00 [0, 0.18]
Mass at nest departure 0.02 [0, 0.21]
Total Conditional 𝑹² 0.67 [0.38, 0.89]
Random effects 𝜎²
Nest 28.21 [5.93, 62.39]
Year 12.22 [0, 68.04]
Residual 20.54 [8.87, 34.02]
Adjusted repeatability 𝑟
Nest 0.49 [0.1, 0.82]
Year 0.20 [0, 0.65]
Residual 0.32 [0.1, 0.68]
Sample sizes 𝑛
Years 6
Individuals 48
Observations (i.e., Nests) 27
Code
mod_flight_age_BC_fits <- 
  as.data.frame(effect(term = "sex", mod = stats_mod_flight_age_BC$mod, 
                       xlevels = list(sex = c("Male", "Female")))) %>%
  mutate(sex_plot = ifelse(sex == "Male", 1.8, 1.2),
         species = "BC")

mod_flight_age_WBC_fits <- 
  as.data.frame(effect(term = "sex", mod = stats_mod_flight_age_WBC$mod, 
                       xlevels = list(sex = c("Male", "Female")))) %>%
  mutate(sex_plot = ifelse(sex == "Male", 1.8, 1.2),
         species = "WBC")

mod_flight_age_fits <- 
  bind_rows(mod_flight_age_BC_fits, 
            mod_flight_age_WBC_fits) %>% 
  `rownames<-`( NULL )

# mod_flight_age_fits[1, 2] - mod_flight_age_fits[2, 2] 
# mod_flight_age_fits[1, 4] - mod_flight_age_fits[2, 4]
# mod_flight_age_fits[1, 5] - mod_flight_age_fits[2, 5] 

flight_age_plot <-
  ggplot2::ggplot() + 
  geom_boxplot(data = flight_dat,
               aes(x = sex_plot, y = flight_age,
                   group = sex, fill = sex),
               color = "grey50",
               width = 0.05, alpha = 0.5,
               position = position_dodge(width = 0)) +
  geom_errorbar(data = mod_flight_age_fits, 
                aes(x = sex_plot, ymax = upper, ymin = lower),
                alpha = 1, color = "black", width = 0.05, lwd = 0.5) + 
  geom_point(data = mod_flight_age_fits, 
             aes(x = sex_plot, y = fit, fill = sex),
             lwd = 3, shape = 21, color= "black") +
  geom_jitter(data = flight_dat, 
              aes(x = sex, y = flight_age, 
                  group = sex, 
                  fill = sex, color = sex), 
              width = 0.02, alpha = 0.2, shape = 19) +
  coord_flip() +
  luke_theme +
  theme(legend.position = "none",
        panel.border = element_blank(),
        strip.background = element_blank(),
        strip.text = element_text(size = 12, face = "italic"),
        axis.title.x = element_text(size = 12),
        # axis.title.y = element_blank(),
        # axis.text.y = element_blank(),
        # axis.ticks.x = element_blank(),
        axis.ticks.y = element_blank(),
        legend.background = element_blank(),
        # panel.grid = element_blank(),
        panel.grid.major.x = element_line(colour = "grey70", size = 0.1),
        axis.title.y = element_blank()) +
  scale_y_continuous(limits = c(4, 42), breaks = seq(5, 40, 5)) +#c(-60, -30, 0, 30, 60)) +
  scale_x_discrete(labels = c("M" = "Male",
                              "F" = "Female")) +
  ylab(expression(paste("Groundling duration (days)" %+-%  "95% CI", sep = ""))) +
  # xlab("Origin") +
  scale_color_manual(values = sex_pal2) +
  scale_fill_manual(values = sex_pal2) +
  facet_grid(species ~ ., labeller = as_labeller(species_names))# +
# annotate(geom = "text", x = 0.5, y = 58,
#          label = "First nests of the season",
#          color = "black", size = 3, fontface = 'italic', hjust = 0)

Code
new_data_WBC <- 
  expand.grid(sex = c("Female","Male"),
              Fledge_age = mean(flight_dat[flight_dat$species == "WBC",]$Fledge_age),
              Fledge_mass = mean(flight_dat[flight_dat$species == "WBC",]$Fledge_mass))

predict_WBC <- 
  predict(stats_mod_flight_age_WBC$mod, 
        newdata = new_data_WBC, 
        re.form = NA, 
        se.fit = TRUE, 
        nsim = 1000)

new_data_BC <- 
  expand.grid(sex = c("Female","Male"),
              Fledge_age = mean(flight_dat[flight_dat$species == "BC",]$Fledge_age),
              Fledge_mass = mean(flight_dat[flight_dat$species == "BC",]$Fledge_mass))

predict_BC <- 
  predict(stats_mod_flight_age_BC$mod, 
        newdata = new_data_BC, 
        re.form = NA, 
        se.fit = TRUE, 
        nsim = 1000)

coucal_flight_age <- 
  data.frame(trait = c("flight_age"),
             species = c("BC", "BC", "WBC", "WBC"),
             sex = c("F", "M", "F", "M"),
             mean = c(predict_BC$fit[1],
                      predict_BC$fit[2],
                      predict_WBC$fit[1],
                      predict_WBC$fit[2]),
             CI_low = c(predict_BC$ci.fit[1, 1],
                        predict_BC$ci.fit[1, 2],
                        predict_WBC$ci.fit[1, 1],
                        predict_WBC$ci.fit[1, 2]),
             CI_high = c(predict_BC$ci.fit[2, 1],
                        predict_BC$ci.fit[2, 2],
                        predict_WBC$ci.fit[2, 1],
                        predict_WBC$ci.fit[2, 2]),
             n_inds = c(filter(flight_dat, species == "BC" & sex == "Female") %>% 
                          summarise(n_ = n_distinct(Ring_ID)) %>% 
                          pull(n_),
                        filter(flight_dat, species == "BC" & sex == "Male") %>% 
                          summarise(n_ = n_distinct(Ring_ID)) %>% 
                          pull(n_),
                        filter(flight_dat, species == "WBC" & sex == "Female") %>% 
                          summarise(n_ = n_distinct(Ring_ID)) %>% 
                          pull(n_),
                        filter(flight_dat, species == "BC" & sex == "Male") %>% 
                          summarise(n_ = n_distinct(Ring_ID)) %>% 
                          pull(n_)),
             n_nests = c(filter(flight_dat, species == "BC" & sex == "Female") %>% 
                           summarise(n_ = n_distinct(nest_ID)) %>% 
                           pull(n_),
                         filter(flight_dat, species == "BC" & sex == "Male") %>% 
                           summarise(n_ = n_distinct(nest_ID)) %>% 
                           pull(n_),
                         filter(flight_dat, species == "WBC" & sex == "Female") %>% 
                           summarise(n_ = n_distinct(nest_ID)) %>% 
                           pull(n_),
                         filter(flight_dat, species == "WBC" & sex == "Male") %>% 
                           summarise(n_ = n_distinct(nest_ID)) %>% 
                           pull(n_)),
             n_years = c(filter(flight_dat, species == "BC" & sex == "Female") %>% 
                           summarise(n_ = n_distinct(year)) %>% 
                           pull(n_),
                         filter(flight_dat, species == "BC" & sex == "Male") %>% 
                           summarise(n_ = n_distinct(year)) %>% 
                           pull(n_),
                         filter(flight_dat, species == "WBC" & sex == "Female") %>% 
                           summarise(n_ = n_distinct(year)) %>% 
                           pull(n_),
                         filter(flight_dat, species == "WBC" & sex == "Male") %>% 
                           summarise(n_ = n_distinct(year)) %>% 
                           pull(n_))) %>% 
  mutate(sd = ifelse(!is.na(CI_low), 
                     approx_sd(x1 = CI_low, x2 = CI_high),
                     CI_low)) 

Bind all parameters together

Code
parameter_distributions <- 
  bind_rows(coucal_egg_survival,
            coucal_clutch_size,
            coucal_HSR,
            coucal_mating_system,
            coucal_fledge_age,
            coucal_flight_age) %>% 
  dplyr::select(trait, species, sex, mean, sd)

Stochastic Demographic Simulation

Code
niter = 10000
set.seed(14)

Black Coucals

pull species- and sex-specific means and sds for each parameter

Code
k_dist_BC = c(pull(filter(coucal_clutch_size, species == "BC"), mean), 
              pull(filter(coucal_clutch_size, species == "BC"), sd))

HSR_dist_BC = c(pull(filter(coucal_HSR, species == "BC"), mean), 
                pull(filter(coucal_HSR, species == "BC"), sd))

h_dist_BC = c(pull(filter(coucal_mating_system, species == "BC"), mean), 
              pull(filter(coucal_mating_system, species == "BC"), sd))

egg_surv_dist_BC = c(pull(filter(coucal_egg_survival, species == "BC"), mean),
                     pull(filter(coucal_egg_survival, species == "BC"), sd))

fledge_age_distF_BC = c(pull(filter(coucal_fledge_age, species == "BC" & sex == "F"), 
                             mean), 
                        pull(filter(coucal_fledge_age, species == "BC" & sex == "F"), 
                             sd))

fledge_age_distM_BC = c(pull(filter(coucal_fledge_age, species == "BC" & sex == "M"), 
                             mean), 
                        pull(filter(coucal_fledge_age, species == "BC" & sex == "M"), 
                             sd))

flight_age_distF_BC = c(pull(filter(coucal_flight_age, species == "BC" & sex == "F"), 
                             mean), 
                        pull(filter(coucal_flight_age, species == "BC" & sex == "F"), 
                             sd))

flight_age_distM_BC = c(pull(filter(coucal_flight_age, species == "BC" & sex == "M"), 
                             mean), 
                        pull(filter(coucal_flight_age, species == "BC" & sex == "M"), 
                             sd))

run stochastic demographic simulation

Code
BC_hazard_ASR_bootstrap_result_w_trans_WBC_ad_surv_stoc <-
  pbsapply(1:niter, run_bootstrap_juv_hazd_ad_surv_ASR,
           offspring = status_dat_all %>% filter(species == "BC"),
           adult_surival_boot_out = trans_WBC_adult_surival_boot_out,
           bootstrap_name = "BC_boot_w_trans_WBC_ad_surv_stoc",
           species = "BC",
           iter_add = 1,
           prefix_number = "boot_w_trans_WBC_ad_surv_stoc_no",
           max_time = 70,
           niter = niter,
           alpha_value = 1.4,
           k_dist = k_dist_BC,
           HSR_dist = HSR_dist_BC,
           h_dist = h_dist_BC,
           egg_surv_dist = egg_surv_dist_BC,
           fledge_age_distF = fledge_age_distF_BC,
           fledge_age_distM = fledge_age_distM_BC,
           flight_age_distF = flight_age_distF_BC,
           flight_age_distM = flight_age_distM_BC, 
           output.directory = "/output/bootstraps/raw")

White-browed Coucals

pull species- and sex-specific means and sds for each parameter

Code
k_dist_WBC = c(pull(filter(coucal_clutch_size, species == "WBC"), mean), 
               pull(filter(coucal_clutch_size, species == "WBC"), sd))

HSR_dist_WBC = c(pull(filter(coucal_HSR, species == "WBC"), mean), 
                 pull(filter(coucal_HSR, species == "WBC"), sd))

h_dist_WBC = c(pull(filter(coucal_mating_system, species == "WBC"), mean), 
               pull(filter(coucal_mating_system, species == "WBC"), sd))

egg_surv_dist_WBC = c(pull(filter(coucal_egg_survival, species == "WBC"), mean),
                      pull(filter(coucal_egg_survival, species == "WBC"), sd))

fledge_age_distF_WBC = c(pull(filter(coucal_fledge_age, species == "WBC" & sex == "F"), 
                              mean), 
                         pull(filter(coucal_fledge_age, species == "WBC" & sex == "F"), 
                              sd))

fledge_age_distM_WBC = c(pull(filter(coucal_fledge_age, species == "WBC" & sex == "M"), 
                              mean), 
                         pull(filter(coucal_fledge_age, species == "WBC" & sex == "M"), 
                              sd))

flight_age_distF_WBC = c(pull(filter(coucal_flight_age, species == "WBC" & sex == "F"), 
                              mean), 
                         pull(filter(coucal_flight_age, species == "WBC" & sex == "F"), 
                              sd))

flight_age_distM_WBC = c(pull(filter(coucal_flight_age, species == "WBC" & sex == "M"),
                              mean), 
                         pull(filter(coucal_flight_age, species == "WBC" & sex == "M"), 
                              sd))

run stochastic demographic simulation

Code
WBC_hazard_ASR_bootstrap_result_w_WBC_ad_surv_stoc <-
  pbsapply(1:niter, run_bootstrap_juv_hazd_ad_surv_ASR,
           offspring = status_dat_all %>% filter(species == "WBC"), 
           adult_surival_boot_out = WBC_adult_surival_boot_out,
           bootstrap_name = "WBC_boot_w_WBC_ad_surv_stoc",
           species = "WBC",
           iter_add = 1,
           prefix_number = "boot_w_WBC_ad_surv_stoc", 
           max_time = 70,
           niter = niter,
           alpha_value = 1.4,
           k_dist = k_dist_WBC,
           HSR_dist = HSR_dist_WBC,
           h_dist = h_dist_WBC,
           egg_surv_dist = egg_surv_dist_WBC,
           fledge_age_distF = fledge_age_distF_WBC,
           fledge_age_distM = fledge_age_distM_WBC,
           flight_age_distF = flight_age_distF_WBC,
           flight_age_distM = flight_age_distM_WBC,
           output.directory = "/output/bootstraps/raw")

Save simulation output

Code
saveRDS(object = BC_hazard_ASR_bootstrap_result_w_trans_WBC_ad_surv_stoc, 
        file = "output/bootstraps/tidy/BC_hazard_ASR_bootstrap_result_w_trans_WBC_ad_surv_stoc.rds")

saveRDS(object = WBC_hazard_ASR_bootstrap_result_w_WBC_ad_surv_stoc, 
        file = "output/bootstraps/tidy/WBC_hazard_ASR_bootstrap_result_w_WBC_ad_surv_stoc.rds")

clean up the output from the bootstrap procedure and save as rds

Code
BC_hazard_rate_boot_tidy <- 
  hazard_boot_out_wrangle(species = "BC", niter = 1000, 
                          output_dir = "output/bootstraps/tidy/",
                          rds_file = 
                           "_hazard_ASR_bootstrap_result_w_trans_WBC_ad_surv_stoc")

WBC_hazard_rate_boot_tidy <- 
  hazard_boot_out_wrangle(species = "WBC", niter = 1000, 
                          output_dir = "output/bootstraps/tidy/",
                          rds_file = 
                            "_hazard_ASR_bootstrap_result_w_WBC_ad_surv_stoc")

save model output

Code
saveRDS(object = BC_hazard_rate_boot_tidy, 
        file = 
          "output/bootstraps/tidy/BC_haz_sur_ASR_boot_tidy_stoc_trans.rds")

saveRDS(object = WBC_hazard_rate_boot_tidy, 
        file = 
          "output/bootstraps/tidy/WBC_haz_sur_ASR_boot_tidy_stoc.rds")

Sex-biased Demographic stages

Visualize the sex-specific variation in cumulative survival and hazard rate of premature individuals derived from the stochastic simulation

Code
grp_means <- 
  bind_rows(BC_hazard_rate_boot_tidy$vital_rate_ests_boot,
          WBC_hazard_rate_boot_tidy$vital_rate_ests_boot) %>% 
  filter(rate == "development") %>% 
  pivot_wider(names_from = stage, values_from = value) %>% 
  mutate(flight_age = fledge_age + flight_age) %>% 
  pivot_longer(c(fledge_age, flight_age), names_to = "stage", values_to = "value") %>% 
  mutate(stage_sex = paste(sex, stage, sep = "_")) %>% 
  group_by(species, stage_sex, sex, stage) %>% 
  dplyr::summarise(grp.mean = mean(value))

BC_density_plotA <- 
  bind_rows(BC_hazard_rate_boot_tidy$vital_rate_ests_boot,
          WBC_hazard_rate_boot_tidy$vital_rate_ests_boot) %>% 
  filter(rate == "development" & species == "BC") %>% 
  pivot_wider(names_from = stage, values_from = value) %>% 
  mutate(flight_age = fledge_age + flight_age) %>% 
  pivot_longer(c(fledge_age, flight_age), names_to = "stage", values_to = "value") %>% 
  mutate(stage_sex = paste(sex, stage, sep = "_")) %>% 
  ggplot() +
  geom_density(aes(value, y=..scaled.., fill = fct_rev(stage_sex)),
               alpha = 0.7, color = "grey40", size = 0.3) + 
  theme_void() +
  theme(legend.position = "none",
        plot.margin = unit(c(0,0,0,0.5), "cm"),
        plot.background = element_rect(fill = 'transparent', color = NA)) +
  scale_fill_manual(values = rev(sex_pal3)) +
  scale_x_continuous(limits = c(0, 70),
                     expand = c(0.01, 0.01)) + 
  ggtitle('A')

BC_density_plotB <- 
  bind_rows(BC_hazard_rate_boot_tidy$vital_rate_ests_boot,
            WBC_hazard_rate_boot_tidy$vital_rate_ests_boot) %>% 
  filter(rate == "development" & species == "BC") %>% 
  pivot_wider(names_from = stage, values_from = value) %>% 
  mutate(flight_age = fledge_age + flight_age) %>% 
  pivot_longer(c(fledge_age, flight_age), names_to = "stage", values_to = "value") %>% 
  mutate(stage_sex = paste(sex, stage, sep = "_")) %>% 
  ggplot() +
  geom_density(aes(value, y=..scaled.., fill = fct_rev(stage_sex)),
               alpha = 0.7, color = "grey40", size = 0.3) + 
  theme_void() +
  theme(legend.position = "none",
        plot.margin = unit(c(0,0,0,0.5), "cm"),
        plot.background = element_rect(fill = 'transparent', color = NA)) +
  scale_fill_manual(values = rev(sex_pal3)) +
  scale_x_continuous(limits = c(0, 70),
                     expand = c(0.01, 0.01)) + 
  ggtitle('B')

WBC_density_plot <-
  bind_rows(BC_hazard_rate_boot_tidy$vital_rate_ests_boot,
            WBC_hazard_rate_boot_tidy$vital_rate_ests_boot) %>% 
  filter(rate == "development" & species == "WBC") %>% 
  pivot_wider(names_from = stage, values_from = value) %>% 
  mutate(flight_age = fledge_age + flight_age) %>% 
  pivot_longer(c(fledge_age, flight_age), names_to = "stage", values_to = "value") %>% 
  mutate(stage_sex = paste(sex, stage, sep = "_")) %>% 
  ggplot() +
  geom_density(aes(value, y=..scaled.., fill = fct_rev(stage_sex)),
               alpha = 0.7, color = "grey40", size = 0.3) + 
  theme_void() +
  theme(legend.position = "none",
        plot.margin = unit(c(0,0,0,0.5), "cm"),
        plot.background = element_rect(fill = 'transparent', color = NA)) +
  scale_fill_manual(values = rev(sex_pal3)) +
  scale_x_continuous(limits = c(0, 70),
                     expand = c(0.01, 0.03))

age_means_BC <- 
  BC_hazard_rate_boot_tidy$hazard_rates_boot %>%
  filter(iter %in% c(as.character(seq(from = 1, to = 1000, by = 1)))) %>% 
  dplyr::group_by(age, sex) %>% 
  dplyr::summarise(mean_rate = mean(estimate, na.rm = TRUE))

age_means_WBC <- 
  WBC_hazard_rate_boot_tidy$hazard_rates_boot %>%
  filter(iter %in% c(as.character(seq(from = 1, to = 1000, by = 1)))) %>% 
  dplyr::group_by(age, sex) %>% 
  dplyr::summarise(mean_rate = mean(estimate, na.rm = TRUE))

BC_surv_plot <-
  ggplot() +
  luke_theme +
  theme(legend.position = "none",
        strip.background = element_blank(),
        strip.text = element_text(size = 12, face = "italic"),
        axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        plot.margin = unit(c(0,0,0.5,0.5), "cm"),
        axis.ticks = element_line(size = 0.25, colour = "grey40"),
        axis.ticks.length = unit(0.1, "cm"),
        axis.text.y  = element_text(size = 8),
        panel.grid.major.x = element_line(colour = "grey70", size = 0.1),
        plot.background = element_rect(fill = 'transparent', color = NA)) +
  geom_vline(data = filter(grp_means, species == "BC" & stage == "fledge_age"),
             aes(xintercept = grp.mean, color = sex), linetype = "dashed", alpha = 1) +
  geom_vline(data = filter(grp_means, species == "BC" & stage == "flight_age"),
             aes(xintercept = grp.mean, color = sex), linetype = "dashed", alpha = 1) +
  geom_line(data = filter(BC_hazard_rate_boot_tidy$hazard_rates_boot, 
                          iter %in% c(as.character(seq(from = 1, to = 1000, by = 1)))),
            aes(x = age, y = estimate, 
                group = interaction(iter, sex), 
                color = sex),
            alpha = 0.05, size = 0.25) +
  geom_line(data = age_means_BC,
            aes(x = age, y = mean_rate, 
                group = sex), 
            color = c(rep("#00496C", 70), rep("#E5BD3A", 70)),
            alpha = 1) +
  facet_grid(species ~ ., labeller = as_labeller(species_names)) +
  scale_colour_manual(values = sex_pal2) +
  ylab("Estimated daily survival rate") +
  xlab("Age (Days since hatching)") + 
  scale_y_continuous(limits = c(0.9, 1),
                     breaks = seq(from = 0.9, to = 1, by = 0.025)) +
  scale_x_continuous(limits = c(0, 70), 
                     breaks = seq(0, 70, by = 5), 
                     expand = c(0.01, 0.01),
                     labels = as.character(seq(0, 70, by = 5))) +
  annotate(geom = "text", y = 0.91, 
           x = mean(c(pull(filter(grp_means, species == "BC" & stage == "fledge_age")[1,5], grp.mean),
                      pull(filter(grp_means, species == "BC" & stage == "fledge_age")[2,5], grp.mean)))/2,
           label = "nestling",
           color = "black", size = 3, fontface = 'italic', hjust = 0.5) +
  annotate(geom = "text", y = 0.91, 
           x = mean(c(pull(filter(grp_means, species == "BC" & stage == "fledge_age")[1,5], grp.mean),
                      pull(filter(grp_means, species == "BC" & stage == "fledge_age")[2,5], grp.mean))) +
             ((mean(c(pull(filter(grp_means, species == "BC" & stage == "flight_age")[1,5], grp.mean),
                      pull(filter(grp_means, species == "BC" & stage == "flight_age")[2,5], grp.mean))) -
                 mean(c(pull(filter(grp_means, species == "BC" & stage == "fledge_age")[1,5], grp.mean),
                        pull(filter(grp_means, species == "BC" & stage == "fledge_age")[2,5], grp.mean))))/2),
           label = "groundling",
           color = "black", size = 3, fontface = 'italic', hjust = 0.5) +
  annotate(geom = "text", y = 0.91, 
           x = mean(c(pull(filter(grp_means, species == "BC" & stage == "flight_age")[1,5], grp.mean),
                      pull(filter(grp_means, species == "BC" & stage == "flight_age")[2,5], grp.mean))) +
             ((70 - mean(c(pull(filter(grp_means, species == "BC" & stage == "flight_age")[1,5], grp.mean), 
                           pull(filter(grp_means, species == "BC" & stage == "flight_age")[2,5], grp.mean))))/2),
           label = "fledgling",
           color = "black", size = 3, fontface = 'italic', hjust = 0.5)

WBC_surv_plot <-
  ggplot() +
  luke_theme +
  theme(legend.position = "none",
        strip.background = element_blank(),
        strip.text = element_text(size = 12, face = "italic"),
        plot.margin = unit(c(0,0,0.5,0.5), "cm"),
        axis.ticks = element_line(size = 0.25, colour = "grey40"),
        axis.ticks.length = unit(0.1, "cm"),
        axis.text.x  = element_text(size = 8),
        axis.text.y  = element_text(size = 8),
        panel.grid.major.x = element_line(colour = "grey70", size = 0.1),
        plot.background = element_rect(fill = 'transparent', color = NA)) +
  geom_vline(data = filter(grp_means, species == "WBC" & stage == "fledge_age"),
             aes(xintercept = grp.mean, color = sex), linetype = "dashed", alpha = 1) +
  geom_vline(data = filter(grp_means, species == "WBC" & stage == "flight_age"),
             aes(xintercept = grp.mean, color = sex), linetype = "dashed", alpha = 1) +
  geom_line(data = filter(WBC_hazard_rate_boot_tidy$hazard_rates_boot, 
                          iter %in% c(as.character(seq(from = 1, to = 1000, by = 1)))),
            aes(x = age, y = estimate, 
                group = interaction(iter, sex), 
                color = sex),
            alpha = 0.05, size = 0.25) +
  geom_line(data = age_means_WBC,
            aes(x = age, y = mean_rate, 
                group = sex), 
                color = c(rep("#00496C", 70), rep("#E5BD3A", 70)),
            alpha = 1) +
  facet_grid(species ~ ., labeller = as_labeller(species_names)) +
  scale_colour_manual(values = sex_pal2) +
  ylab("                                        Estimated daily survival rate") +
  xlab("Age (days since hatching)") + 
  scale_y_continuous(limits = c(0.9, 1),
                     breaks = seq(from = 0.9, to = 1, by = 0.025)) +
  scale_x_continuous(limits = c(0, 70), 
                     breaks = seq(0, 70, by = 5), 
                     expand = c(0.01, 0.01),
                     labels = as.character(seq(0, 70, by = 5))) +
  annotate(geom = "text", y = 0.91, 
           x = mean(c(pull(filter(grp_means, species == "WBC" & stage == "fledge_age")[1,5], grp.mean),
                      pull(filter(grp_means, species == "WBC" & stage == "fledge_age")[2,5], grp.mean)))/2,
           label = "nestling",
           color = "black", size = 3, fontface = 'italic', hjust = 0.5) +
  annotate(geom = "text", y = 0.91, 
           x = mean(c(pull(filter(grp_means, species == "WBC" & stage == "fledge_age")[1,5], grp.mean),
                      pull(filter(grp_means, species == "WBC" & stage == "fledge_age")[2,5], grp.mean))) +
             ((mean(c(pull(filter(grp_means, species == "WBC" & stage == "flight_age")[1,5], grp.mean),
                      pull(filter(grp_means, species == "WBC" & stage == "flight_age")[2,5], grp.mean))) -
                 mean(c(pull(filter(grp_means, species == "WBC" & stage == "fledge_age")[1,5], grp.mean),
                        pull(filter(grp_means, species == "WBC" & stage == "fledge_age")[2,5], grp.mean))))/2),
           label = "groundling",
           color = "black", size = 3, fontface = 'italic', hjust = 0.5) +
  annotate(geom = "text", y = 0.91, 
           x = mean(c(pull(filter(grp_means, species == "WBC" & stage == "flight_age")[1,5], grp.mean),
                      pull(filter(grp_means, species == "WBC" & stage == "flight_age")[2,5], grp.mean))) +
             ((70 - mean(c(pull(filter(grp_means, species == "WBC" & stage == "flight_age")[1,5], grp.mean), 
                           pull(filter(grp_means, species == "WBC" & stage == "flight_age")[2,5], grp.mean))))/2),
           label = "fledgling",
           color = "black", size = 3, fontface = 'italic', hjust = 0.5)

# BC_coxme <- coxme(Surv(exit, event) ~ sex + (1|year) + (1|nest_ID), data =  status_dat_all %>% filter(species == "BC"))  
# summary(BC_coxme)

BC_survfit_sex <- survfit(Surv(exit, event) ~ sex , data = status_dat_all %>% filter(species == "BC"))

BC_plot <- 
  ggsurvplot(BC_survfit_sex,
             risk.table = F,
             ncensor.plot = F,
             conf.int = T,
             ylab = "Cumulative Survival",
             xlab = "Age (Days since hatching)",
             xlim = c(0,70),
             break.time.by = 10,
             palette = sex_pal2,
             ggtheme = luke_theme + theme(plot.title = element_text(hjust = 0.5, face = "italic"),
                                          legend.background = element_rect(fill="transparent")),
             legend =  c(0.7, 0.9),
             font.legend = c(10, "plain", "black"),
             legend.title = "",
             legend.labs = c("Female", "Male"))

# WBC_coxme <- coxme(Surv(exit, event) ~ sex + (1|year) + (1|nest_ID), data =  status_dat_all %>% filter(species == "WBC"))  
# summary(WBC_coxme)

WBC_survfit_sex <- survfit(Surv(exit, event) ~ sex , data = status_dat_all %>% filter(species == "WBC"))

WBC_plot <- 
  ggsurvplot(WBC_survfit_sex,
             risk.table = F,
             ncensor.plot = F,
             conf.int = T,
             ylab = "Cumulative Survival",
             xlab = "Age (Days since hatching)",
             xlim = c(0,70),
             break.time.by = 10,
             palette = sex_pal2,
             ggtheme = luke_theme + theme(plot.title = element_text(hjust = 0.5, face = "italic"),
                                          legend.background = element_rect(fill="transparent")),
             legend =  c(0.7, 0.9),
             font.legend = c(10, "plain", "black"),
             legend.title = "",
             legend.labs = c("Female", "Male"))

BC_plot1 <- 
  BC_plot$plot +
  geom_vline(data = filter(grp_means, species == "BC" & stage == "fledge_age"),
             aes(xintercept = grp.mean, color = sex), linetype = "dashed", alpha = 1, color = c("#be9c2e", "#016392")) +
  geom_vline(data = filter(grp_means, species == "BC" & stage == "flight_age"),
             aes(xintercept = grp.mean, color = sex), linetype = "dashed", alpha = 1, color = c("#be9c2e", "#016392")) +
  luke_theme +
  theme(legend.position = c(0.8, 0.8),
        plot.margin = unit(c(0,0,0.5,0), "cm"),
        legend.text=element_text(size = 10),
        legend.title = element_blank(),
        legend.key.size = unit(0.3, 'cm'),
        # panel.border = element_blank(),
        strip.background = element_blank(),
        strip.text = element_text(size = 12, face = "italic"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        # axis.ticks.y = element_blank(),
        axis.ticks = element_line(size = 0.25, colour = "grey40"),
        axis.ticks.length = unit(0.1, "cm"),
        axis.text.x  = element_blank(),
        axis.text.y  = element_text(size = 8),
        panel.grid.major.x = element_line(colour = "grey70", size = 0.1),
        # plot.title = element_text(hjust = 0.5, face = "italic"),
        legend.background = element_rect(fill="transparent"),
        plot.background = element_rect(fill = 'transparent', color = NA)) +
  scale_x_continuous(limits = c(0, 70), 
                     breaks = seq(0, 70, by = 5), 
                     expand = c(0.01, 0.01),
                     labels = as.character(seq(0, 70, by = 5))) +#c(-60, -30, 0, 30, 60)) +
  # ylab(expression(paste("Cumulative survival" %+-%  "95% CI", sep = ""))) +
  # xlab("Age (Days since hatching)") +
  scale_color_manual(values = sex_pal2,
                     labels = c("Female", "Male")) +
  scale_fill_manual(values = sex_pal2,
                    labels = c("Female", "Male")) +
  annotate(geom = "text", y = 0.08, 
           x = mean(c(pull(filter(grp_means, species == "WBC" & stage == "fledge_age")[1,5], grp.mean),
                      pull(filter(grp_means, species == "WBC" & stage == "fledge_age")[2,5], grp.mean)))/2,
           label = "nestling",
           color = "black", size = 3, fontface = 'italic', hjust = 0.5) +
  annotate(geom = "text", y = 0.08, 
           x = mean(c(pull(filter(grp_means, species == "WBC" & stage == "fledge_age")[1,5], grp.mean),
                      pull(filter(grp_means, species == "WBC" & stage == "fledge_age")[2,5], grp.mean))) +
             ((mean(c(pull(filter(grp_means, species == "WBC" & stage == "flight_age")[1,5], grp.mean),
                      pull(filter(grp_means, species == "WBC" & stage == "flight_age")[2,5], grp.mean))) -
                 mean(c(pull(filter(grp_means, species == "WBC" & stage == "fledge_age")[1,5], grp.mean),
                        pull(filter(grp_means, species == "WBC" & stage == "fledge_age")[2,5], grp.mean))))/2),
           label = "groundling",
           color = "black", size = 3, fontface = 'italic', hjust = 0.5) +
  annotate(geom = "text", y = 0.08, 
           x = mean(c(pull(filter(grp_means, species == "WBC" & stage == "flight_age")[1,5], grp.mean),
                      pull(filter(grp_means, species == "WBC" & stage == "flight_age")[2,5], grp.mean))) +
             ((70 - mean(c(pull(filter(grp_means, species == "WBC" & stage == "flight_age")[1,5], grp.mean), 
                           pull(filter(grp_means, species == "WBC" & stage == "flight_age")[2,5], grp.mean))))/2),
           label = "fledgling",
           color = "black", size = 3, fontface = 'italic', hjust = 0.5)

WBC_plot1 <- 
  WBC_plot$plot +
  geom_vline(data = filter(grp_means, species == "WBC" & stage == "fledge_age"),
             aes(xintercept = grp.mean, color = sex), linetype = "dashed", alpha = 1, color = c("#be9c2e", "#016392")) +
  geom_vline(data = filter(grp_means, species == "WBC" & stage == "flight_age"),
             aes(xintercept = grp.mean, color = sex), linetype = "dashed", alpha = 1, color = c("#be9c2e", "#016392")) +
  luke_theme +
  theme(legend.position = "none",
        plot.margin = unit(c(0,0,0.5,0), "cm"),
        # panel.border = element_blank(),
        strip.background = element_blank(),
        strip.text = element_text(size = 12, face = "italic"),
        axis.title.x = element_text(size = 12),
        axis.ticks = element_line(size = 0.25, colour = "grey40"),
        axis.ticks.length = unit(0.1, "cm"),
        axis.text.x  = element_text(size = 8),
        axis.text.y  = element_text(size = 8),
        # axis.ticks.y = element_blank(),
        legend.background = element_blank(),
        panel.grid.major.x = element_line(colour = "grey70", size = 0.1),
        plot.background = element_rect(fill = 'transparent', color = NA)) +
        # axis.title.y = element_blank(),
        # axis.text.y = element_blank()) +
  scale_x_continuous(limits = c(0, 70), 
                     breaks = seq(0, 70, by = 5), 
                     expand = c(0.01, 0.01),
                     labels = as.character(seq(0, 70, by = 5))) +#c(-60, -30, 0, 30, 60)) +
  # xlab("Age (Days since hatching)") +
  scale_color_manual(values = sex_pal2) +
  scale_fill_manual(values = sex_pal2) +
  ylab("                                    Cumulative survival ± 95% CI") +
  xlab("Age (days since hatching)") +
  annotate(geom = "text", y = 0.08, 
           x = mean(c(pull(filter(grp_means, species == "WBC" & stage == "fledge_age")[1,5], grp.mean),
                      pull(filter(grp_means, species == "WBC" & stage == "fledge_age")[2,5], grp.mean)))/2,
           label = "nestling",
           color = "black", size = 3, fontface = 'italic', hjust = 0.5) +
  annotate(geom = "text", y = 0.08, 
           x = mean(c(pull(filter(grp_means, species == "WBC" & stage == "fledge_age")[1,5], grp.mean),
                      pull(filter(grp_means, species == "WBC" & stage == "fledge_age")[2,5], grp.mean))) +
             ((mean(c(pull(filter(grp_means, species == "WBC" & stage == "flight_age")[1,5], grp.mean),
                      pull(filter(grp_means, species == "WBC" & stage == "flight_age")[2,5], grp.mean))) -
                 mean(c(pull(filter(grp_means, species == "WBC" & stage == "fledge_age")[1,5], grp.mean),
                        pull(filter(grp_means, species == "WBC" & stage == "fledge_age")[2,5], grp.mean))))/2),
           label = "groundling",
           color = "black", size = 3, fontface = 'italic', hjust = 0.5) +
  annotate(geom = "text", y = 0.08, 
           x = mean(c(pull(filter(grp_means, species == "WBC" & stage == "flight_age")[1,5], grp.mean),
                      pull(filter(grp_means, species == "WBC" & stage == "flight_age")[2,5], grp.mean))) +
             ((70 - mean(c(pull(filter(grp_means, species == "WBC" & stage == "flight_age")[1,5], grp.mean), 
                           pull(filter(grp_means, species == "WBC" & stage == "flight_age")[2,5], grp.mean))))/2),
           label = "fledgling",
           color = "black", size = 3, fontface = 'italic', hjust = 0.5)

Cox-proportional hazards model

Here we statistically test for sex-specific age-dependent mortality hazard, while accounting for non-independence among siblings by clustering by nest_ID. We also test the assumptions of proportional hazards based on scaled Schoenfeld residuals (met: black coucal: P = 0.5; white-browed coucal: P = 0.08). In summary, we find no sex effect on hazard rate in either species (black coucal: eβ = 0.91, 95% CI = 0.70–1.19, P = 0.49; white-browed coucal: eβ = 0.90, 95% CI = 0.67–1.19, P = 0.47).

Code
BC_cox_ph_sex <- 
  coxph(Surv(time = entry, 
             time2 = exit, 
             event = event) ~ sex + cluster(nest_ID), 
        data = status_dat_all %>% filter(species == "BC"))

WBC_cox_ph_sex <- 
  coxph(Surv(time = entry, 
             time2 = exit, 
             event = event) ~ sex + cluster(nest_ID), 
        data = status_dat_all %>% filter(species == "WBC"))

Black Coucal offspring

no sex effect on hazard rate: eβ = 0.91, 95% CI = 0.70–1.19, P = 0.493

Code
summary(BC_cox_ph_sex)
Call:
coxph(formula = Surv(time = entry, time2 = exit, event = event) ~ 
    sex, data = status_dat_all %>% filter(species == "BC"), cluster = nest_ID)

  n= 636, number of events= 208 

         coef exp(coef) se(coef) robust se      z Pr(>|z|)
sexM -0.09308   0.91112  0.13911   0.13589 -0.685    0.493

     exp(coef) exp(-coef) lower .95 upper .95
sexM    0.9111      1.098    0.6981     1.189

Concordance= 0.503  (se = 0.019 )
Likelihood ratio test= 0.45  on 1 df,   p=0.5
Wald test            = 0.47  on 1 df,   p=0.5
Score (logrank) test = 0.45  on 1 df,   p=0.5,   Robust = 0.47  p=0.5

  (Note: the likelihood ratio and score tests assume independence of
     observations within a cluster, the Wald and robust score tests do not).

Schoenfeld residuals and diagnostic plot show that assumptions are met (P = 0.5)

Code
cox.zph(BC_cox_ph_sex)
       chisq df   p
sex    0.448  1 0.5
GLOBAL 0.448  1 0.5
Code
plot(cox.zph(BC_cox_ph_sex, transform = "identity"))

White-browed Coucal offspring

no sex effect on hazard rate: eβ = 0.90, 95% CI = 0.67–1.19, P = 0.446

Code
summary(WBC_cox_ph_sex)
Call:
coxph(formula = Surv(time = entry, time2 = exit, event = event) ~ 
    sex, data = status_dat_all %>% filter(species == "WBC"), 
    cluster = nest_ID)

  n= 435, number of events= 167 

        coef exp(coef) se(coef) robust se      z Pr(>|z|)
sexM -0.1105    0.8954   0.1549    0.1451 -0.762    0.446

     exp(coef) exp(-coef) lower .95 upper .95
sexM    0.8954      1.117    0.6737      1.19

Concordance= 0.532  (se = 0.021 )
Likelihood ratio test= 0.51  on 1 df,   p=0.5
Wald test            = 0.58  on 1 df,   p=0.4
Score (logrank) test = 0.51  on 1 df,   p=0.5,   Robust = 0.58  p=0.4

  (Note: the likelihood ratio and score tests assume independence of
     observations within a cluster, the Wald and robust score tests do not).

Schoenfeld residuals and diagnostic plot show that assumptions are met (P = 0.081)

Code
cox.zph(WBC_cox_ph_sex)
       chisq df     p
sex     3.04  1 0.081
GLOBAL  3.04  1 0.081
Code
plot(cox.zph(WBC_cox_ph_sex, transform = "identity"))

Visualising offspring survival dynamics

Survival dynamics of radio-marked Black (top) and White-browed (bottom) Coucal offspring over the first 70 days since hatching. A) Kaplan-Meier plots of the sex-specific cumulative survival. B) Inverted sex-specific smoothed hazard functions (right) showing the estimated daily survival rate. Fine lines show the estimated function of survival for a bootstrap that randomly sampled one offspring per brood to account for non-independence. Thick lines show the average trend for the bootstrapping procedure. Density plots show the sex-specific posterior distributions of age-points of nest departure and first flight given the bootstrap simulation.

Visualize the sex- and stage-specific variation in the survival rates derived from the stochastic simulation

Code
# calculate the sex differences in stage specific rates
BC_sex_diff_hazard_output <- 
  sex_diff_hazard(boot_out_list = BC_hazard_rate_boot_tidy, 
                  niter = 1000) %>% 
  mutate(species = "BC")

WBC_sex_diff_hazard_output <- 
  sex_diff_hazard(boot_out_list = WBC_hazard_rate_boot_tidy, 
                  niter = 1000) %>% 
  mutate(species = "WBC")

# consolidate results
sex_diff_survival_output <- 
  dplyr::bind_rows(BC_sex_diff_hazard_output,
                   WBC_sex_diff_hazard_output) %>% 
  dplyr::filter(stage %in% c("HSR", "Nestling", "Groundling", 
                             "Fledgling", "Adult")) %>% 
  mutate(stage = factor(stage, 
                        levels = c("HSR", "Nestling", "Groundling", 
                                   "Fledgling", "Adult")),
         species = factor(species, 
                          levels = c("BC", "WBC")),
         difference = -difference)

# check that there are no nosensical values
# sex_diff_survival_output %>% 
#   filter(difference < -1 & difference > 1)
# 
# ggplot(data = filter(sex_diff_survival_output, stage != "h")) +
#   geom_histogram(aes(difference), binwidth = 0.01) +
#   facet_grid(stage ~ species) +
#   scale_x_continuous(limits = c(-0.5, 0.5)) +
#   geom_vline(xintercept = 0)

# calculate some summary statistics
sex_diff_survival_summary <- 
  sex_diff_survival_output %>%
  dplyr::group_by(stage, species) %>%
  dplyr::summarise(avg = mean(difference, na.rm = TRUE),
                   median = median(difference, na.rm = TRUE),
                   var = var(difference, na.rm = TRUE),
                   max = max(difference, na.rm = TRUE),
                   min = min(difference, na.rm = TRUE),
                   sd = sd(difference, na.rm = TRUE))

junk0 <- 
  data.frame(species = c("BC", "WBC"),
             stage = "HSR",
             difference = c(NA,NA))

sex_diff_survival_output2 <- 
  bind_rows(sex_diff_survival_output, junk0) %>% 
  mutate(difference = as.numeric(difference),
         stage = factor(stage, 
                        levels = c("HSR", "Nestling", "Groundling", 
                                   "Fledgling", "Adult")),
         species = factor(species, 
                          levels = c("BC", "WBC")))

# Figure 2a: plot the sex-biases in survival across the three stages
vital_rate_theme <- 
  theme_bw() +
  theme(
    text = element_text(family = "Verdana"),
    legend.position = "none",
    panel.background = element_rect(fill = "transparent", colour = NA),
    plot.background = element_rect(fill = "transparent", colour = NA),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.ticks.length = unit(0.1, "cm"),
    panel.border = element_blank(),
    panel.spacing = unit(0.3, "lines"),
    strip.background = element_blank()
  )

surv_diff_plot <-
  ggplot(aes(y = difference, x = stage, fill = stage), 
         data = slice_sample(group_by(sex_diff_survival_output, species, stage), n = 1000)) +
  geom_violin(draw_quantiles = c(0.025, 0.5, 0.975), color = "grey10",
              scale = "width", trim = TRUE, adjust = 1, size = 0.25) +
  vital_rate_theme +
  facet_grid(. ~ species, labeller = as_labeller(species_names)) +
  theme(axis.title.x = element_text(size = 12, colour = "black"),
        axis.text.x  = element_text(size = 8, angle = 45, hjust = 1, vjust = 1, 
                                    colour = "black"),
        axis.ticks.x = element_line(size = 0.2),
        axis.title.y = element_text(size = 12, hjust = 0.5, vjust = -3),
        axis.text.y  = element_text(size = 8, colour = "black"),
        axis.ticks.y = element_line(size = 0.2),
        plot.margin = unit(c(0.2, 1.35, 0.42, 0.3), "cm"),
        strip.text = element_text(size = 12, colour = "black", face = "italic")
  ) +
  scale_fill_manual(values = cbPalette_sex_diff) +
  scale_y_continuous(limits = c(-0.65, 0.65),
                     breaks = c(-0.5, -0.25, 0, 0.25, 0.5),
                     expand = c(0, 0), position = "left",
                     labels = c("-0.50","-0.25",
                                expression(phantom("-")*"0.00"),
                                expression(phantom("-")*"0.25"),
                                expression(phantom("-")*"0.50"))) +
  xlab("Life stage") +
  ylab("Sex bias\n") +
  scale_x_discrete(labels = c(
    "HSR" = "hatching\nsex ratio",
    "Nestling" = "nestling\nsurvival",
    "Groundling" = "groundling\nsurvival",
    "Fledgling" = "fledgling\nsurvival",
    "Adult" = "adult\nsurvival"))

background <-
  ggplot(aes(y = difference, x = stage, fill = stage),
         data = sex_diff_survival_output) +
  annotate("rect", xmin = -Inf, xmax = Inf, 
           ymin = -Inf, ymax = 0, alpha = 0.7,
           fill = pull(ggthemes_data$wsj$palettes$colors6[3,2])) +
  annotate("rect", xmin = -Inf, xmax = Inf, 
           ymin = 0, ymax = Inf, alpha = 0.7,
           fill = pull(ggthemes_data$wsj$palettes$colors6[2,2])) +
  annotate("text", x = c(5), y = c(-0.5),
           label = c("\u2640"), size = 10,
           family = "Menlo",
           vjust = c(0), hjust = c(0.5), colour = "grey10") +
  annotate("text", x = c(1), y = c(0.5),
           label = c("\u2642"), size = 10,
           family = "Menlo",
           vjust = c(1), hjust = c(0.5), colour = "grey10") +
  facet_grid(. ~ species, labeller = as_labeller(species_names)) +
  vital_rate_theme +
  theme(axis.title.x = element_text(size = 12, colour = "white"),
        axis.text.x  = element_text(size = 8, angle = 45, hjust = 1, 
                                    vjust = 1, colour = "white"),
        axis.ticks.x = element_line(size = 0.2, colour = "white"),
        axis.title.y = element_text(size = 12, hjust = 0.5, vjust = -1, 
                                    colour = "white"),
        axis.text.y  = element_text(size = 8, colour = "white"),
        axis.ticks.y = element_line(size = 0.2, colour = "white"),
        plot.margin = unit(c(0.2, 1.35, 2.05, 0.3), "cm"),
        strip.text = element_text(size = 12, colour = "white")
  ) +
  scale_x_discrete(labels = c(
    "HSR" = "hatching\nsex ratio",
    "Nestling" = "nestling\nsurvival",
    "Groundling" = "groundling\nsurvival",
    "Fledgling" = "fledgling\nsurvival",
    "Adult" = "adult\nsurvival")) +
  scale_y_continuous(limits = c(-0.65, 0.65), 
                     breaks = c(-0.5, -0.25, 0, 0.25, 0.5), 
                     expand = c(0, 0), position = "left",
                     labels = c("-0.50","-0.25", 
                                expression(phantom("-")*"0.00"), 
                                expression(phantom("-")*"0.25"),
                                expression(phantom("-")*"0.50"))) +
  xlab("Life stage") +
  ylab("Sex bias\n")

Derived Adult Sex Ratio

Code
ASR_boot <- 
  bind_rows(BC_hazard_rate_boot_tidy$ASR_ests_boot,
            WBC_hazard_rate_boot_tidy$ASR_ests_boot) %>%
  mutate(species = factor(species, levels = c("BC", "WBC")))

CI <- 0.95

ASR_boot_summary <- 
  ASR_boot %>%
  dplyr::group_by(species) %>%
  dplyr::summarise(ucl_ASR = stats::quantile(M_Adult, (1 - CI)/2, na.rm = TRUE),
                   lcl_ASR = stats::quantile(M_Adult, 1 - (1 - CI)/2, na.rm = TRUE),
                   avg_ASR = mean(M_Adult),
                   med_ASR = median(M_Adult),
                   max_ASR = max(M_Adult),
                   min_ASR = min(M_Adult))

Figure_2b <- 
  ggplot() +
  annotate("rect", xmin=0, xmax=0.5, ymin=0, ymax=610, alpha=0.6,
           fill= pull(ggthemes_data$wsj$palettes$colors6[3,2])) +
  annotate("rect", xmin=0.5, xmax=1, ymin=0, ymax=610, alpha=0.6,
           fill= pull(ggthemes_data$wsj$palettes$colors6[2,2])) +
  annotate("text", x = c(0.05), y = c(305),
           label = c("\u2640"), size = 10, colour = "grey10",
           family="Menlo", vjust = c(0), hjust = c(0.5)) +
  annotate("text", x = c(0.95), y = c(305),
           label = c("\u2642"), size = 10, colour = "grey10",
           family="Menlo", vjust = c(1), hjust = c(0.5)) +
  geom_histogram(binwidth = 0.01, data = ASR_boot, aes(x = M_Adult), fill = "grey30") +
  geom_errorbarh(data = ASR_boot_summary, aes(y = 600, x = lcl_ASR, xmin = lcl_ASR, xmax = ucl_ASR), 
                 color = "black", size = 0.3, linetype = "solid") +
  coord_flip() +
  facet_grid(. ~ species, labeller = as_labeller(species_names)) +
  theme_bw() +
  theme(text = element_text(family="Verdana"),
        legend.position = "none",
        axis.title.x = element_text(size=12, vjust=-0.1),
        axis.text.x  = element_text(size=8, angle = 45, hjust = 1, colour = "black"),
        axis.title.y = element_text(size=12, hjust=0.5, vjust = -2.75, margin = margin(0, 11, 0, 0)),
        axis.text.y  = element_text(size=8, colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.y = element_line(size = 0.2, colour = "grey40"),
        axis.ticks.length = unit(0.1, "cm"),
        axis.ticks.x = element_line(size = 0.2, colour = "grey40"),
        panel.border = element_blank(),
        plot.margin = unit(c(0.2, 1.35, 0.405, 0.08), "cm"),
        panel.spacing = unit(0.3, "lines"),
        strip.background = element_blank(),
        strip.text = element_text(size = 12, colour = "black", face = "italic")) +
  ylab("Bootstrap frequency") +
  xlab("Adult sex ratio\n(proportion male, 95% CI)") +
  scale_x_continuous(limits = c(0, 1), expand = c(0, 0)) +
  scale_y_continuous(limits = c(0, 610), expand = c(0, 0), breaks=c(0, 100, 200, 300, 400, 500))

Life-table Response Experiment

Functions for the LTRE analysis

sensitivity_analysis()

Sensitivity analysis function takes the vital rate summary of the bootstrap procedure and conducts a perturbation analysis on each rate to assess how a proportional change in a given vital rate changes the ASR

arguments:

  • vital_rate_summary: object (a list of all vital rates) produced by the make_mprime_matrix() and make_treat_matrix() functions below
  • matrix_str: the matrix_structure expression
  • h: the mating system index specified within vital_rate_summary
  • k: the clutch size specified within vital_rate_summary
  • HSR: the hatching sex ratio specified within vital_rate_summary
  • niter: number of time steps to project matrix model
  • ASR: the stable stage distribution ASR derived from the matrix_ASR() function applied to the mprime or treat matrix
  • lambda: the stable stage distribution lambda derived from the matrix_ASR() function applied to the mprime or treat matrix
Code
sensitivity_analysis <-
  function(vital_rate_summary, matrix_str, h = 1, k = 4, 
           HSR, 
           niter = 1000, ASR, lambda){
    
    # make a list of all parameters
    vr <-
      list(F_Nestling_survival = vital_rate_summary$F_Nestling_survival,
           F_Groundling_survival = vital_rate_summary$F_Groundling_survival,
           F_Fledgling_survival = vital_rate_summary$F_Fledgling_survival,
           F_Adult_survival = vital_rate_summary$F_Adult_survival,
           M_Nestling_survival = vital_rate_summary$M_Nestling_survival,
           M_Groundling_survival = vital_rate_summary$M_Groundling_survival,
           M_Fledgling_survival = vital_rate_summary$M_Fledgling_survival,
           M_Adult_survival = vital_rate_summary$M_Adult_survival)
    
    # number of stages in the matrix
    no_stages <- sqrt(length(matrix_str))
    
    # Define plover life-stages of the Ceuta snowy plover matrix model
    stages <- c("F_1st_year",  "F_Adult",  "M_1st_year",  "M_Adult")
    
    # an empty t by x matrix
    stage <- matrix(numeric(no_stages * niter), nrow = no_stages)
    
    # an empty t vector to store the population sizes
    pop <- numeric(niter)
    
    # dataframe to store the perturbation results
    ASR_pert_results <-
      data.frame(parameter = c("F_Nestling_survival", "F_Groundling_survival", 
                               "F_Fledgling_survival", "F_Adult_survival",
                               "M_Nestling_survival", "F_Groundling_survival", 
                               "M_Fledgling_survival", "M_Adult_survival",
                               "h", "k", "HSR", "ISR"),
                 sensitivities = numeric(length(vr) + 4),
                 elasticities = numeric(length(vr) + 4))
    
    lambda_pert_results <-
      data.frame(parameter = c("F_Nestling_survival", "F_Groundling_survival", 
                               "F_Fledgling_survival", "F_Adult_survival",
                               "M_Nestling_survival", "F_Groundling_survival", 
                               "M_Fledgling_survival", "M_Adult_survival",
                               "h", "k", "HSR", "ISR"),
                 sensitivities = numeric(length(vr) + 4),
                 elasticities = numeric(length(vr) + 4))
    
    # specifiy how many survival rates there are
    n <- length(vr)
    
    # create vectors of perturbations to test on parameters of the matrix model
    vr_nums <- seq(0, 1, 0.01) # proportional changes in survival and HSR (i.e., between 0 an 1)
    h_nums <- seq(0, 2, 0.02) # proportional changes in h index (i.e., between 0 and 2)
    k_nums <- seq(3, 5, 0.02) # proportional changes in k (i.e, between 3 and 5)
    
    # create empty dataframes to store the perturbation results for ASR and lambda
    vr_pert_ASR <- matrix(numeric(n * length(vr_nums)),
                          ncol = n, dimnames = list(vr_nums, names(vr)))
    
    h_pert_ASR <- matrix(numeric(length(h_nums)),
                         ncol = 1, dimnames = list(h_nums, "h"))
    
    k_pert_ASR <- matrix(numeric(length(k_nums)),
                         ncol = 1, dimnames = list(k_nums, "k"))
    
    HSR_pert_ASR <- matrix(numeric(length(vr_nums)),
                           ncol = 1, dimnames = list(vr_nums, "HSR"))
    
    vr_pert_lambda <- matrix(numeric(n * length(vr_nums)),
                             ncol = n, dimnames = list(vr_nums, names(vr)))
    
    h_pert_lambda <- matrix(numeric(length(h_nums)),
                            ncol = 1, dimnames = list(h_nums, "h"))
    
    k_pert_lambda <- matrix(numeric(length(k_nums)),
                            ncol = 1, dimnames = list(k_nums, "k"))
    
    HSR_pert_lambda <- matrix(numeric(length(vr_nums)),
                              ncol = 1, dimnames = list(vr_nums, "HSR"))
    
    #### perturbation of survival rates ####
    for (g in 1:n) { # pick a column (i.e., a variable) 
      vr2 <- vr # reset the vital rates to the original
      
      for (i in 1:length(vr_nums)) # pick a perturbation level
      {
        
        vr2[[g]] <- vr_nums[i] # specify the vital rate with the new perturbation level
        
        A <- matrix(sapply(matrix_str, eval, vr2, NULL), 
                    nrow = sqrt(length(matrix_str)), byrow=TRUE, 
                    dimnames = list(stages, stages)) # build the matrix with the new value
        
        # reset the starting stage distribution for simulation (all with 10 individuals)
        m <- rep(10, no_stages) 
        
        for (j in 1:niter) { # project the matrix through t iteration
          # stage distribution at time t
          stage[,j] <- m
          
          # population size at time t
          pop[j] <- sum(m)
          
          # number of male adults at time t
          M2 <- stage[4, j] * A["M_Adult", "M_Adult"]
          
          # number of female adults at time t
          F2 <- stage[2, j] * A["F_Adult", "F_Adult"]
          
          # Female freq-dep fecundity of Female chicks
          A["F_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h)) * (1 - HSR) )
          
          # Female freq-dep fecundity of Male chicks
          A["M_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h)) * HSR)
          
          # Male freq-dep fecundity of Female chicks
          A["F_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h)) * (1 - HSR))
          
          # Male freq-dep fecundity of Male chicks
          A["M_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h)) * HSR)
          
          # define the new n (i.e., new stage distribution at time t)
          m <- A %*% (m)
        }
        
        # define rownames of stage matrix
        rownames(stage) <- rownames(A)
        
        # define colnames of stage matrix
        colnames(stage) <- 0:(niter - 1)
        
        # calculate the proportional stable stage distribution
        stage <- apply(stage, 2, function(x) x / sum(x))
        
        # define stable stage as the last stage
        stable.stage <- stage[, niter]
        
        # calc ASR as the proportion of the adult stable stage class that is male
        vr_pert_ASR[i, g] <- stable.stage[no_stages] / (stable.stage[no_stages/2] + 
                                                          stable.stage[no_stages])
        
        # calc lambda as the pop change in the counts of the last two iterations
        vr_pert_lambda[i, g] <- pop[niter]/pop[niter - 1]
      }
      
      # get the spline function of ASR
      spl_ASR <- 
        smooth.spline(vr_pert_ASR[which(!is.na(vr_pert_ASR[, g])), g] ~ 
                        names(vr_pert_ASR[which(!is.na(vr_pert_ASR[, g])), g]))
      
      # estimate the slope of the tangent of the spline at the vital rate
      ASR_pert_results[g, 2] <- predict(spl_ASR, x = vr[[g]], deriv = 1)$y
      
      # re-scale sensitivity into elasticity
      ASR_pert_results[g, 3] <- vr[[g]] / ASR * ASR_pert_results[g, 2]
      
      # do the same steps but for lambda
      spl_lambda <- 
        smooth.spline(vr_pert_lambda[which(!is.na(vr_pert_lambda[, g])), g] ~ 
                        names(vr_pert_lambda[which(!is.na(vr_pert_lambda[, g])), g]))
      
      lambda_pert_results[g, 2] <- predict(spl_lambda, x = vr[[g]], deriv = 1)$y
      
      lambda_pert_results[g, 3] <- vr[[g]] / lambda * lambda_pert_results[g, 2]
    }
    
    #### perturbation of the h index parameter ####
    for (i in 1:length(h_nums)) # pick a perturbation level
    {
      A <- matrix(sapply(matrix_str, eval, vr, NULL), 
                  nrow = sqrt(length(matrix_str)), byrow=TRUE, 
                  dimnames = list(stages, stages)) 
      
      # reset the starting stage distribution for simulation (all with 10 individuals)
      m <- rep(10, no_stages) 
      
      for (j in 1:niter) { # project the matrix through t iteration
        # stage distribution at time t
        stage[,j] <- m
        
        # population size at time t
        pop[j] <- sum(m)
        
        # number of male adults at time t
        M2 <- stage[4, j] * A["M_Adult", "M_Adult"]
        
        # number of female adults at time t
        F2 <- stage[2, j] * A["F_Adult", "F_Adult"]
        
        # Female freq-dep fecundity of Female chicks
        A["F_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h_nums[i])) * (1 - HSR) )
        
        # Female freq-dep fecundity of Male chicks
        A["M_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h_nums[i])) * HSR)
        
        # Male freq-dep fecundity of Female chicks
        A["F_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h_nums[i])) * (1 - HSR))
        
        # Male freq-dep fecundity of Male chicks
        A["M_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h_nums[i])) * HSR)
        
        # define the new n (i.e., new stage distribution at time t)
        m <- A %*% (m)
      }
      # define rownames of stage matrix
      rownames(stage) <- rownames(A)
      
      # define colnames of stage matrix
      colnames(stage) <- 0:(niter - 1)
      
      # calculate the proportional stable stage distribution
      stage <- apply(stage, 2, function(x) x / sum(x))
      
      # define stable stage as the last stage
      stable.stage <- stage[, niter]
      
      # calc ASR as the proportion of the adult stable stage class that is male
      h_pert_ASR[i,] <- stable.stage[no_stages] / (stable.stage[no_stages / 2] + 
                                                     stable.stage[no_stages])
      
      # calc lambda as the pop change in the counts of the last two iterations
      h_pert_lambda[i, ] <- pop[niter]/pop[niter - 1]
      
    }
    # get the spline function of ASR
    spl_ASR <- 
      smooth.spline(h_pert_ASR[which(!is.na(h_pert_ASR)), 1] ~ 
                      names(h_pert_ASR[which(!is.na(h_pert_ASR)), ]))
    
    # estimate the slope of the tangent of the spline at the vital rate
    ASR_pert_results[n + 1, 2] <- predict(spl_ASR, x = h, deriv = 1)$y
    
    # re-scale sensitivity into elasticity
    ASR_pert_results[n + 1, 3] <- h / ASR * ASR_pert_results[n + 1, 2]
    
    # do the same steps but for lambda
    spl_lambda <- 
      smooth.spline(h_pert_lambda[which(!is.na(h_pert_lambda)), 1] ~ 
                      names(h_pert_lambda[which(!is.na(h_pert_lambda)), ]))
    lambda_pert_results[n + 1, 2] <- predict(spl_lambda, x = h, deriv = 1)$y
    lambda_pert_results[n + 1, 3] <- h / lambda * lambda_pert_results[n + 1, 2]
    
    #### perturbation of k parameter ####
    for (i in 1:length(k_nums)) # pick a perturbation level
    {
      A <- matrix(sapply(matrix_str, eval, vr, NULL), 
                  nrow = sqrt(length(matrix_str)), byrow=TRUE, 
                  dimnames = list(stages, stages))
      
      # reset the starting stage distribution for simulation (all with 10 individuals)
      m <- rep(10, no_stages) 
      
      for (j in 1:niter) { # project the matrix through t iteration
        # stage distribution at time t
        stage[,j] <- m
        
        # population size at time t
        pop[j] <- sum(m)
        
        # number of male adults at time t
        M2 <- stage[4, j] * A["M_Adult", "M_Adult"]
        
        # number of female adults at time t
        F2 <- stage[2, j] * A["F_Adult", "F_Adult"]
        
        # Female freq-dep fecundity of Female chicks
        A["F_1st_year", "F_Adult"] <- ((k_nums[i] * M2) / (M2 + (F2 / h)) * (1 - HSR) )
        
        # Female freq-dep fecundity of Male chicks
        A["M_1st_year", "F_Adult"] <- ((k_nums[i] * M2) / (M2 + (F2 / h)) * HSR)
        
        # Male freq-dep fecundity of Female chicks
        A["F_1st_year", "M_Adult"] <- ((k_nums[i] * F2) / (M2 + (F2 / h)) * (1 - HSR))
        
        # Male freq-dep fecundity of Male chicks
        A["M_1st_year", "M_Adult"] <- ((k_nums[i] * F2) / (M2 + (F2 / h)) * HSR)
        
        # define the new n (i.e., new stage distribution at time t)
        m <- A %*% (m)
      }
      # define rownames of stage matrix
      rownames(stage) <- rownames(A)
      
      # define colnames of stage matrix
      colnames(stage) <- 0:(niter - 1)
      
      # calculate the proportional stable stage distribution
      stage <- apply(stage, 2, function(x) x / sum(x))
      
      # define stable stage as the last stage
      stable.stage <- stage[, niter]
      
      # calc ASR as the proportion of the adult stable stage class that is male
      k_pert_ASR[i,] <- stable.stage[no_stages] / (stable.stage[no_stages/2] + 
                                                     stable.stage[no_stages])
      
      # calc lambda as the pop change in the counts of the last two iterations
      k_pert_lambda[i, ] <- pop[niter]/pop[niter - 1]
    }
    
    # get the spline function of ASR
    spl_ASR <- 
      smooth.spline(k_pert_ASR[which(!is.na(k_pert_ASR)), 1] ~ 
                      names(k_pert_ASR[which(!is.na(k_pert_ASR)), ]))
    
    # estimate the slope of the tangent of the spline at the vital rate
    
    ASR_pert_results[n + 2, 2] <- predict(spl_ASR, x = k, deriv = 1)$y
    
    # re-scale sensitivity into elasticity
    ASR_pert_results[n + 2, 3] <- k / ASR * ASR_pert_results[n + 2, 2]
    
    # do the same steps but for lambda
    spl_lambda <- 
      smooth.spline(h_pert_lambda[which(!is.na(h_pert_lambda)), 1] ~ 
                      names(h_pert_lambda[which(!is.na(h_pert_lambda)), ]))
    lambda_pert_results[n + 2, 2] <- predict(spl_lambda, x = k, deriv = 1)$y
    lambda_pert_results[n + 2, 3] <- k / lambda * lambda_pert_results[n + 2, 2]
    
    #### perturbation of HSR ####
    for (i in 1:length(vr_nums)) # pick a perturbation level
    {
      A <- matrix(sapply(matrix_str, eval, vr, NULL), 
                  nrow = sqrt(length(matrix_str)), byrow=TRUE, 
                  dimnames = list(stages, stages))
      
      # reset the starting stage distribution for simulation (all with 10 individuals)
      m <- rep(10, no_stages) 
      
      for (j in 1:niter) { # project the matrix through t iteration
        # stage distribution at time t
        stage[,j] <- m
        
        # population size at time t
        pop[j] <- sum(m)
        
        # number of male adults at time t
        M2 <- stage[4, j] * A["M_Adult", "M_Adult"]
        
        # number of female adults at time t
        F2 <- stage[2, j] * A["F_Adult", "F_Adult"]
        
        # Female freq-dep fecundity of Female chicks
        A["F_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h)) * (1 - vr_nums[i]) )
        
        # Female freq-dep fecundity of Male chicks
        A["M_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h)) * vr_nums[i])
        
        # Male freq-dep fecundity of Female chicks
        A["F_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h)) * (1 - vr_nums[i]))
        
        # Male freq-dep fecundity of Male chicks
        A["M_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h)) * vr_nums[i])
        
        # define the new n (i.e., new stage distribution at time t)
        m <- A %*% (m)
        
      }
      # define rownames of stage matrix
      rownames(stage) <- rownames(A)
      
      # define colnames of stage matrix
      colnames(stage) <- 0:(niter - 1)
      
      # calculate the proportional stable stage distribution
      stage <- apply(stage, 2, function(x) x / sum(x))
      
      # define stable stage as the last stage
      stable.stage <- stage[, niter]
      
      # calc ASR as the proportion of the adult stable stage class that is male
      HSR_pert_ASR[i,] <- stable.stage[no_stages] / (stable.stage[no_stages/2] + 
                                                       stable.stage[no_stages])
      
      # calc lambda as the pop change in the counts of the last two iterations
      HSR_pert_lambda[i, ] <- pop[niter] / pop[niter - 1]
      
    }
    # get the spline function of ASR
    spl_ASR <- 
      smooth.spline(HSR_pert_ASR[which(!is.na(HSR_pert_ASR)), 1] ~ 
                      names(HSR_pert_ASR[which(!is.na(HSR_pert_ASR)), ]))
    
    # estimate the slope of the tangent of the spline at the vital rate    
    ASR_pert_results[n + 3, 2] <- predict(spl_ASR, x = HSR, deriv = 1)$y
    
    # re-scale sensitivity into elasticity
    ASR_pert_results[n + 3, 3] <- HSR / ASR * ASR_pert_results[n + 3, 2]
    
    # do the same steps but for lambda
    spl_lambda <- 
      smooth.spline(h_pert_lambda[which(!is.na(h_pert_lambda)), 1] ~ 
                      names(h_pert_lambda[which(!is.na(h_pert_lambda)), ]))
    lambda_pert_results[n + 3, 2] <- predict(spl_lambda, x = HSR, deriv = 1)$y
    lambda_pert_results[n + 3, 3] <- HSR/lambda * lambda_pert_results[n + 3, 2]
    
    #### store all results into a list ----
    result <- list(ASR_pert_results = ASR_pert_results,
                   lambda_pert_results = lambda_pert_results,
                   HSR_pert_ASR = HSR_pert_ASR,
                   HSR_pert_lambda = HSR_pert_lambda,
                   k_pert_ASR = k_pert_ASR,
                   k_pert_lambda = k_pert_lambda,
                   h_pert_ASR = h_pert_ASR,
                   h_pert_lambda = h_pert_lambda,
                   vr_pert_ASR = vr_pert_ASR,
                   vr_pert_lambda = vr_pert_lambda)
    
  }

matrix_structure

matrix_structure is an expression that determines the structure of the two-sex matrix from a set of vital rates

Code
matrix_structure <- expression(
  # top row of matrix
  0, NA, 0, NA,
  
  # second row of matrix
  (F_Nestling_survival * F_Groundling_survival * F_Fledgling_survival), F_Adult_survival, 0, 0,
  
  # third row of matrix
  0, NA, 0, NA,
  
  # fourth row of matrix
  0, 0, (M_Nestling_survival * M_Groundling_survival * M_Fledgling_survival), M_Adult_survival
  
)

make_treat_matrix()

This function makes a treatment matrix from a summary of the bootstrapped vital rates

arguments:

  • survival_rates_boot_summary: the $vital_rate_ests_boot element in the tidied _hazard_rate_boot_tidy object from above (a summary of the bootstrapped vital rates
  • species: species name (either “BC” or “WBC”)
  • h: the mating system index specified in the parameter_distributions
  • k: the clutch size specified in the parameter_distributions
  • HSR: the hatching sex ratio specified in the parameter_distributions
Code
make_treat_matrix <- 
  function(survival_rates_boot_summary, 
           species_name, 
           h, 
           k, 
           HSR){
    
    list(F_Nestling_survival = filter(survival_rates_boot_summary,
                                      species == species_name & vital_rate == "Female_nestling_survival")[, 4],
         F_Groundling_survival = filter(survival_rates_boot_summary,
                                        species == species_name & vital_rate == "Female_groundling_survival")[, 4],
         F_Fledgling_survival = filter(survival_rates_boot_summary,
                                       species == species_name & vital_rate == "Female_fledgling_survival")[, 4],
         F_Adult_survival = filter(survival_rates_boot_summary,
                                   species == species_name & vital_rate == "Female_adult_survival")[, 4],
         M_Nestling_survival = filter(survival_rates_boot_summary,
                                      species == species_name & vital_rate == "Male_nestling_survival")[, 4],
         M_Groundling_survival = filter(survival_rates_boot_summary,
                                        species == species_name & vital_rate == "Male_groundling_survival")[, 4],
         M_Fledgling_survival = filter(survival_rates_boot_summary,
                                       species == species_name & vital_rate == "Male_fledgling_survival")[, 4],
         M_Adult_survival = filter(survival_rates_boot_summary,
                                   species == species_name & vital_rate == "Male_adult_survival")[, 4],
         Egg_survival = filter(survival_rates_boot_summary,
                               species == species_name & vital_rate == "NA_egg_survival")[, 4],
         
         # Define h (harem size, h = 1 is monogamy) and k (clutch size)
         h = h,
         k = k,
         
         # Define primary sex ratio
         HSR = HSR)
  }

make_mprime_matrix()

This function makes a prime-matrix (i.e., a matrix halfway between the treatment matrix and the unbiased matrix) from a set of vital rate summaries

arguments:

  • survival_rates_boot_summary: the $vital_rate_ests_boot element in the tidied _hazard_rate_boot_tidy object from above (a summary of the bootstrapped vital rates
  • species: species name (either “BC” or “WBC”)
  • h: the mating system index specified in the parameter_distributions
  • k: the clutch size specified in the parameter_distributions
  • HSR: the hatching sex ratio specified in the parameter_distributions
  • sex: the sex that you wish to specify as the baseline vital rates for the LTRE
Code
make_mprime_matrix <- 
  function(survival_rates_boot_summary, 
           species_name, 
           h, 
           k, 
           HSR, 
           sex){
    if(sex == "male"){
      list(F_Nestling_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Female_nestling_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Male_nestling_survival")[, 4]) / 2,
           F_Groundling_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Female_groundling_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Male_groundling_survival")[, 4]) / 2,
           F_Fledgling_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Female_fledgling_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Male_fledgling_survival")[, 4]) / 2,
           F_Adult_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Female_adult_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Male_adult_survival")[, 4]) / 2,
           M_Nestling_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Male_nestling_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Male_nestling_survival")[, 4]) / 2,
           M_Groundling_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Male_groundling_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Male_groundling_survival")[, 4]) / 2,
           M_Fledgling_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Male_fledgling_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Male_fledgling_survival")[, 4]) / 2,
           M_Adult_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Male_adult_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Male_adult_survival")[, 4]) / 2,
           Egg_survival = filter(survival_rates_boot_summary,
                                 species == species_name)[15, 4],
           
           # Define h (harem size, h = 1 is monogamy) and k (clutch size)
           h = (h + 1) / 2,
           k = k,
      
           # Define primary sex ratio
           HSR = (HSR + 0.5) / 2)
           

    }
    else{
      list(F_Nestling_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Female_nestling_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Female_nestling_survival")[, 4]) / 2,
           F_Groundling_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Female_groundling_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Female_groundling_survival")[, 4]) / 2,
           F_Fledgling_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Female_fledgling_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Female_fledgling_survival")[, 4]) / 2,
           F_Adult_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Female_adult_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Female_adult_survival")[, 4]) / 2,
           M_Nestling_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Female_nestling_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Male_nestling_survival")[, 4]) / 2,
           M_Groundling_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Female_groundling_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Male_groundling_survival")[, 4]) / 2,
           M_Fledgling_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Female_fledgling_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Male_fledgling_survival")[, 4]) / 2,
           M_Adult_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Female_adult_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Male_adult_survival")[, 4]) / 2,
           Egg_survival = filter(survival_rates_boot_summary,
                                 species == species_name)[15, 4],
           
           # Define h (harem size, h = 1 is monogamy) and k (clutch size)
           h = (h + 1) / 2,
           k = k,
           
           # Define primary sex ratio
           HSR = (HSR + 0.5) / 2)
    }
  }

LTRE_analysis()

LTRE_analysis() compares the relative difference in ASR response of an unbiased two-sex matrix vs. the observed sex-specific matrix

arguments:

  • Mprime_sens: the object produced by the sensitivity_analysis() applied to the Mprime matrix
  • matrix_str: the object of the matrix_structure expression
  • vital_rates: the object produced by the make_treat_matrix()
  • species_name: species name (either “BC” or “WBC”)
  • sex: the sex that you wish to specify as the baseline vital rates for the LTRE
Code
LTRE_analysis <-
  function(Mprime_sens, matrix_str, vital_rates, species_name, sex){
    
    # make empty dataframes to store LTRE results for ASR and lambda
    LTRE_ASR <-
      data.frame(parameter = c("Nestling survival", "Groundling survival",
                               "Fledgling survival", "Adult survival",
                               "Hatching sex ratio", 
                               "Mating system"),
                 contribution = numeric(6))
    LTRE_lambda <-
      data.frame(parameter = c("Nestling survival", "Groundling survival",
                               "Fledgling survival", "Adult survival",
                               "Hatching sex ratio",
                               "Mating system"),
                 contribution = numeric(6))
    
    # run a for loop to extract the parameter contributions
    if (sex == "male"){
      # male rates scenario
      for(i in 1:nrow(LTRE_ASR))
      {
        LTRE_ASR[i, 2] <-
          
          # survival rates
          ifelse(i < 5, (vital_rates[[i + 4]] - vital_rates[[i]]) * 
                   Mprime_sens$ASR_pert_results$sensitivities[i + 4],
                 
                 # HSR
                 ifelse(i == 5, (vital_rates[[12]] - 0.5) * 
                          Mprime_sens$ASR_pert_results$sensitivities[11],
                               
                               # mating system
                               (1 - vital_rates[[10]]) * 
                                 Mprime_sens$ASR_pert_results$sensitivities[9]))
      }
      for(i in 1:nrow(LTRE_lambda))
      {
        # survival rates
        ifelse(i < 5, (vital_rates[[i + 4]] - vital_rates[[i]]) * 
                 Mprime_sens$lambda_pert_results$sensitivities[i + 4],
               
               # HSR
               ifelse(i == 5, (vital_rates[[12]] - 0.5) * 
                        Mprime_sens$lambda_pert_results$sensitivities[11],
                             
                             # mating system
                             (1 - vital_rates[[10]]) * 
                               Mprime_sens$lambda_pert_results$sensitivities[9]))
      }
    }
    else{
      # female rates scenario
      for(i in 1:nrow(LTRE_ASR))
      {
        LTRE_ASR[i, 2] <-
          
          # survival rates
          ifelse(i < 5, (vital_rates[[i]] - vital_rates[[i + 4]]) * 
                   Mprime_sens$ASR_pert_results$sensitivities[i],
                 
                 # HSR
                 ifelse(i == 5, (vital_rates[[12]] - 0.5) * 
                          Mprime_sens$ASR_pert_results$sensitivities[11],
                        
                               # mating system
                               (1 - vital_rates[[10]]) * 
                                 Mprime_sens$ASR_pert_results$sensitivities[9]))
      }
      
      for(i in 1:nrow(LTRE_lambda))
      {
        # survival rates
        ifelse(i < 5, (vital_rates[[i]] - vital_rates[[i + 4]]) * 
                 Mprime_sens$lambda_pert_results$sensitivities[i],
               
               # HSR
               ifelse(i == 5, (vital_rates[[12]] - 0.5) * 
                        Mprime_sens$lambda_pert_results$sensitivities[11],
                             
                             # mating system
                             (1 - vital_rates[[10]]) * 
                               Mprime_sens$lambda_pert_results$sensitivities[9]))
      }
    }
    LTRE_ASR$parameter <- 
      factor(LTRE_ASR$parameter, levels = c("Adult survival",
                                            "Fledgling survival",
                                            "Groundling survival",
                                            "Nestling survival",
                                            "Hatching sex ratio",
                                            "Mating system"))
    LTRE_lambda$parameter <- 
      factor(LTRE_lambda$parameter, levels = c("Adult survival",
                                               "Fledgling survival",
                                               "Groundling survival",
                                               "Nestling survival",
                                               "Hatching sex ratio",
                                               "Mating system"))
    
    LTRE_ASR$model <- "ASR"
    LTRE_ASR$species <- species_name
    
    LTRE_lambda$model <- "lambda"
    LTRE_lambda$species <- species_name
    
    LTRE_results <- list(LTRE_ASR = LTRE_ASR,
                         LTRE_lambda = LTRE_lambda)
    
    LTRE_results
  }

LTRE Analysis

Here we conduct the Life Table Response Experiment to decompose the difference in response between two or more “treatments” by weighting the difference in parameter values by the parameter’s contribution to the response (i.e., its sensitivity), and summing over all parameters. Our LTRE compares the observed scenario (M), to a null scenario with no sex differences (M0), whereby all male survival rates were set equal to the female rates (M0♀), the hatching sex ratio was unbiased (i.e., ρ = 0.5), and mating system was monogamous (i.e., h = 1). Thus, our LTRE identifies the drivers of ASR bias by decomposing the absolute parameter contributions to the difference between the ASR predicted by our model and an unbiased ASR.

Code
BC_hazard_rate_boot_tidy$vital_rate_ests_boot$iter <- 
  as.factor(BC_hazard_rate_boot_tidy$vital_rate_ests_boot$iter)
WBC_hazard_rate_boot_tidy$vital_rate_ests_boot$iter <- 
  as.factor(WBC_hazard_rate_boot_tidy$vital_rate_ests_boot$iter)

# summarize the vital rates
survival_rates_boot_summary <-
  bind_rows(BC_hazard_rate_boot_tidy$vital_rate_ests_boot, WBC_hazard_rate_boot_tidy$vital_rate_ests_boot) %>% 
  mutate(vital_rate = paste(sex, stage, rate, sep = "_")) %>% 
  Rmisc::summarySE(.,
                   measurevar = "value",
                   groupvars = c("vital_rate", "species"),
                   conf.interval = 0.95) %>% 
  arrange(species, vital_rate)


BC_VR_treat <- 
  make_treat_matrix(survival_rates_boot_summary = survival_rates_boot_summary,
                    h = 1/pull(filter(parameter_distributions, species == "BC" & trait == "mating_system"), mean),
                    HSR = pull(filter(parameter_distributions, species == "BC" & trait == "hatching_sex_ratio"), mean),
                    k = pull(filter(parameter_distributions, species == "BC" & trait == "clutch_size"), mean),
                    species_name = "BC")

BC_VR_mprime_male <- 
  make_mprime_matrix(survival_rates_boot_summary = survival_rates_boot_summary,
                     h = 1/pull(filter(parameter_distributions, species == "BC" & trait == "mating_system"), mean),
                     HSR = pull(filter(parameter_distributions, species == "BC" & trait == "hatching_sex_ratio"), mean),
                     k = pull(filter(parameter_distributions, species == "BC" & trait == "clutch_size"), mean),
                     species_name = "BC",
                     sex = "male")

BC_VR_mprime_female <- 
  make_mprime_matrix(survival_rates_boot_summary = survival_rates_boot_summary,
                     h = 1/pull(filter(parameter_distributions, species == "BC" & trait == "mating_system"), mean),
                     HSR = pull(filter(parameter_distributions, species == "BC" & trait == "hatching_sex_ratio"), mean),
                     k = pull(filter(parameter_distributions, species == "BC" & trait == "clutch_size"), mean),
                     species_name = "BC",
                     sex = "female")
WBC_VR_treat <- 
  make_treat_matrix(survival_rates_boot_summary = survival_rates_boot_summary,
                    h = 1/pull(filter(parameter_distributions, species == "WBC" & trait == "mating_system"), mean),
                    HSR = pull(filter(parameter_distributions, species == "WBC" & trait == "hatching_sex_ratio"), mean),
                    k = pull(filter(parameter_distributions, species == "WBC" & trait == "clutch_size"), mean),
                    species_name = "WBC")

WBC_VR_mprime_male <- 
  make_mprime_matrix(survival_rates_boot_summary = survival_rates_boot_summary,
                     h = 1/pull(filter(parameter_distributions, species == "WBC" & trait == "mating_system"), mean),
                     HSR = pull(filter(parameter_distributions, species == "WBC" & trait == "hatching_sex_ratio"), mean),
                     k = pull(filter(parameter_distributions, species == "WBC" & trait == "clutch_size"), mean),
                     species_name = "WBC",
                     sex = "male")

WBC_VR_mprime_female <- 
  make_mprime_matrix(survival_rates_boot_summary = survival_rates_boot_summary,
                     h = 1/pull(filter(parameter_distributions, species == "WBC" & trait == "mating_system"), mean),
                     HSR = pull(filter(parameter_distributions, species == "WBC" & trait == "hatching_sex_ratio"), mean),
                     k = pull(filter(parameter_distributions, species == "WBC" & trait == "clutch_size"), mean),
                     species_name = "WBC",
                     sex = "female")

BC_treatment_matrix <- coucal_matrix(BC_VR_treat)
BC_M_prime_matrix_male <- coucal_matrix(BC_VR_mprime_male)
BC_M_prime_matrix_female <- coucal_matrix(BC_VR_mprime_female)

WBC_treatment_matrix <- coucal_matrix(WBC_VR_treat)
WBC_M_prime_matrix_male <- coucal_matrix(WBC_VR_mprime_male)
WBC_M_prime_matrix_female <- coucal_matrix(WBC_VR_mprime_female)

BC_treatment_ASR_analysis <- 
  matrix_ASR(M = BC_treatment_matrix,  
             h = BC_VR_treat$h, 
             k = BC_VR_treat$k,
             HSR = BC_VR_treat$HSR,
             iterations = 100, 
             num_boot = 1, 
             iter_add = 1,
             species = "BC")

BC_ASR_treat <- BC_treatment_ASR_analysis$ASR

WBC_treatment_ASR_analysis <- 
  matrix_ASR(M = WBC_treatment_matrix,  
             h = WBC_VR_treat$h, 
             k = WBC_VR_treat$k,
             HSR = WBC_VR_treat$HSR,
             iterations = 100, 
             num_boot = 1, 
             iter_add = 1,
             species = "BC")

WBC_ASR_treat <- WBC_treatment_ASR_analysis$ASR

BC_M_prime_ASR_analysis_male <- 
  matrix_ASR(M = BC_M_prime_matrix_male, 
             h = BC_VR_mprime_male$h, 
             k = BC_VR_mprime_male$k,
             HSR = BC_VR_mprime_male$HSR,
             iterations = 100,
             num_boot = 1,
             iter_add = 1,
             species = "BC")

BC_ASR_mprime_male <- BC_M_prime_ASR_analysis_male$ASR

BC_M_prime_ASR_analysis_female <- 
  matrix_ASR(M = BC_M_prime_matrix_female, 
             h = BC_VR_mprime_female$h, 
             k = BC_VR_mprime_female$k,
             HSR = BC_VR_mprime_female$HSR,
             iterations = 100,
             num_boot = 1,
             iter_add = 1,
             species = "BC")

BC_ASR_mprime_female <- BC_M_prime_ASR_analysis_female$ASR

WBC_M_prime_ASR_analysis_male <- 
  matrix_ASR(M = WBC_M_prime_matrix_male, 
             h = WBC_VR_mprime_male$h, 
             k = WBC_VR_mprime_male$k,
             HSR = WBC_VR_mprime_male$HSR,
             iterations = 100,
             num_boot = 1,
             iter_add = 1,
             species = "BC")

WBC_ASR_mprime_male <- WBC_M_prime_ASR_analysis_male$ASR

WBC_M_prime_ASR_analysis_female <- 
  matrix_ASR(M = WBC_M_prime_matrix_female, 
             h = WBC_VR_mprime_female$h, 
             k = WBC_VR_mprime_female$k,
             HSR = WBC_VR_mprime_female$HSR,
             iterations = 100,
             num_boot = 1,
             iter_add = 1,
             species = "BC")

WBC_ASR_mprime_female <- WBC_M_prime_ASR_analysis_female$ASR

BC_lambda_treat <- BC_treatment_ASR_analysis$lambda

BC_lambda_mprime_male <- BC_M_prime_ASR_analysis_male$lambda

BC_lambda_mprime_female <- BC_M_prime_ASR_analysis_female$lambda

WBC_lambda_treat <- WBC_treatment_ASR_analysis$lambda

WBC_lambda_mprime_male <- WBC_M_prime_ASR_analysis_male$lambda

WBC_lambda_mprime_female <- WBC_M_prime_ASR_analysis_female$lambda

BC_treat_sensitivity_analysis <- 
  sensitivity_analysis(vital_rate_summary = BC_VR_treat, 
                       matrix_str = matrix_structure, 
                       h = BC_VR_treat$h, 
                       k = BC_VR_treat$k, 
                       HSR = BC_VR_treat$HSR, 
                       niter = 100, 
                       ASR = BC_ASR_treat,
                       lambda = BC_lambda_treat)

BC_Mprime_sensitivity_analysis_male <- 
  sensitivity_analysis(vital_rate_summary = BC_VR_mprime_male, 
                       matrix_str = matrix_structure, 
                       h = BC_VR_mprime_male$h, 
                       k = BC_VR_mprime_male$k, 
                       HSR = BC_VR_mprime_male$HSR, 
                       niter = 100, 
                       ASR = BC_ASR_mprime_male,
                       lambda = BC_lambda_mprime_male)

BC_Mprime_sensitivity_analysis_female <- 
  sensitivity_analysis(vital_rate_summary = BC_VR_mprime_female, 
                       matrix_str = matrix_structure, 
                       h = BC_VR_mprime_female$h, 
                       k = BC_VR_mprime_female$k, 
                       HSR = BC_VR_mprime_female$HSR, 
                       niter = 100, 
                       ASR = BC_ASR_mprime_female,
                       lambda = BC_lambda_mprime_female) 

WBC_treat_sensitivity_analysis <- 
  sensitivity_analysis(vital_rate_summary = WBC_VR_treat, 
                       matrix_str = matrix_structure, 
                       h = WBC_VR_treat$h, 
                       k = WBC_VR_treat$k, 
                       HSR = WBC_VR_treat$HSR, 
                       niter = 100, 
                       ASR = WBC_ASR_treat,
                       lambda = WBC_lambda_treat)

WBC_Mprime_sensitivity_analysis_male <- 
  sensitivity_analysis(vital_rate_summary = WBC_VR_mprime_male, 
                       matrix_str = matrix_structure, 
                       h = WBC_VR_mprime_male$h, 
                       k = WBC_VR_mprime_male$k, 
                       HSR = WBC_VR_mprime_male$HSR, 
                       niter = 100, 
                       ASR = WBC_ASR_mprime_male,
                       lambda = WBC_lambda_mprime_male) 

WBC_Mprime_sensitivity_analysis_female <- 
  sensitivity_analysis(vital_rate_summary = WBC_VR_mprime_female, 
                       matrix_str = matrix_structure, 
                       h = WBC_VR_mprime_female$h, 
                       k = WBC_VR_mprime_female$k, 
                       HSR = WBC_VR_mprime_female$HSR, 
                       niter = 100, 
                       ASR = WBC_ASR_mprime_female,
                       lambda = WBC_lambda_mprime_female)
Code
BC_LTRE_male <- 
  LTRE_analysis(Mprime_sens = BC_Mprime_sensitivity_analysis_male, 
                matrix_str = matrix_str, 
                vital_rates = BC_VR_treat,
                species_name = "BC",
                sex = "male")

BC_LTRE_female <- 
  LTRE_analysis(Mprime_sens = BC_Mprime_sensitivity_analysis_female, 
                matrix_str = matrix_str, 
                vital_rates = BC_VR_treat,
                species_name = "BC",
                sex = "female")

WBC_LTRE_male <- 
  LTRE_analysis(Mprime_sens = WBC_Mprime_sensitivity_analysis_male, 
                matrix_str = matrix_str, 
                vital_rates = WBC_VR_treat,
                species_name = "WBC",
                sex = "male")

WBC_LTRE_female <- 
  LTRE_analysis(Mprime_sens = WBC_Mprime_sensitivity_analysis_female, 
                matrix_str = matrix_str, 
                vital_rates = WBC_VR_treat,
                species_name = "WBC",
                sex = "female")

LTRE_coucal_male_ASR <- 
  bind_rows(BC_LTRE_male$LTRE_ASR, WBC_LTRE_male$LTRE_ASR) %>% 
  mutate(sex = "male")

LTRE_coucal_female_ASR <- 
  bind_rows(BC_LTRE_female$LTRE_ASR, WBC_LTRE_female$LTRE_ASR) %>% 
  mutate(sex = "female")

LTRE_coucal_ASR <- rbind(LTRE_coucal_male_ASR, LTRE_coucal_female_ASR)

LTRE_coucal_ASR$parameter <- 
  factor(LTRE_coucal_ASR$parameter, 
         levels = c("Hatching sex ratio",
                    "Nestling survival",
                    "Groundling survival",
                    "Fledgling survival",
                    "Adult survival",
                    "Mating system"))

LTRE Results

Visualization of the LTRE results

Code
LTRE_ASR_plot_background <-
  ggplot2::ggplot(data = LTRE_coucal_ASR,
                  aes(x = parameter, y = contribution, fill = parameter)) +
  annotate("rect", xmin=-Inf, xmax=Inf, ymin=-0.15, ymax=0, alpha=0.6,
           fill= sex_pal2[1]) + 
  annotate("rect", xmin=-Inf, xmax=Inf, ymin=0, ymax=0.15, alpha=0.6,
           fill= sex_pal2[2]) + 
  annotate("text", x = c(2), y = c(-0.14),
           label = c("\u2640"), size = 10, family = "Menlo",
           vjust = c(0), hjust = c(0.5), colour = "grey10") +
  annotate("text", x = c(2), y = c(0.14),
           label = c("\u2642"), size = 10, family = "Menlo",
           vjust = c(1), hjust = c(0.5), colour = "grey10") +
  facet_grid(sex ~ species, 
             labeller = labeller(.cols = species_names, .rows = analysis_names)) +
  vital_rate_theme +
  theme(axis.title.x = element_text(size = 12, vjust = -0.1, colour = "white"),
        axis.text.x  = element_text(size = 8, angle = 45, hjust = 1, colour = "white"),
        axis.ticks.x = element_line(size = 0.2, colour = "white"),
        strip.text.x = element_text(size = 12, colour = "white"),
        axis.title.y = element_text(size = 12, hjust=0.5, vjust = 3.5, colour = "white"),
        axis.text.y  = element_text(size = 8, colour = "white"),
        axis.ticks.y = element_blank(),
        strip.text.y = element_text(size = 12, colour = "white"),
        plot.margin = unit(c(0.2, 0.85,
                             0.78, 0.6), "cm")) +
  
  scale_x_discrete(labels = c("Hatching sex ratio" = expression(italic("\u03C1")),
                              "Nestling survival" = expression(italic(S["n"])),
                              "Groundling survival" = expression(italic(S["g"])),
                              "Fledgling survival" = expression(italic(S["f"])),
                              "Adult survival" = expression(italic(phi["a"])),
                              "Mating system" = expression(italic("h")))) +
  scale_y_continuous(limits = c(-0.15 , 0.15), expand = c(0, 0),
                     breaks = c(-0.15, -0.1, -0.05, 0, 0.05, 0.1, 0.15)) +
  ylab("Absolute contribution to adult sex ratio bias") +
  xlab("Sex bias in parameter")

LTRE_ASR_plot <-
  ggplot2::ggplot() +
  geom_bar(data = LTRE_coucal_ASR,
           aes(x = parameter, y = contribution, fill = parameter), 
           color = "black", stat = "identity", size = 0.2) + 
  facet_grid(sex ~ species, 
             labeller = labeller(.cols = species_names, .rows = analysis_names)) +
  vital_rate_theme +
  theme(axis.title.x = element_text(size=12, vjust=-0.1),
        axis.text.x  = element_text(size=8, angle = 45, hjust = 1),
        axis.ticks.x = element_line(size = 0.2, colour = "grey40"),
        strip.text = element_text(size = 12, colour = "black", face = "italic"),
        panel.background = element_rect(fill = "transparent", colour = NA),
        plot.background = element_rect(fill = "transparent", colour = NA),
        axis.title.y = element_text(size=12, hjust=0.5, vjust = 3.5),
        axis.text.y  = element_text(size=8),
        axis.ticks.y = element_line(size = 0.2, colour = "grey40"),
        strip.text.y = element_text(size = 12),
        plot.margin = unit(c(0.2, 0.85, 0.2, 0.6), "cm")) +
  scale_fill_manual(values = cbPalette_LTRE) +
  scale_y_continuous(limits = c(-0.15 ,0.15), expand = c(0, 0),
                     breaks = c(-0.15, -0.1, -0.05, 0, 0.05, 0.1, 0.15)) +
  scale_x_discrete(labels = c("Hatching sex ratio" = expression(italic("\u03C1")),
                              "Nestling survival" = expression(italic(S["n"])),
                              "Groundling survival" = expression(italic(S["g"])),
                              "Fledgling survival" = expression(italic(S["f"])),
                              "Adult survival" = expression(italic(phi["a"])),
                              "Mating system" = expression(italic("h")))) +
  ylab("Absolute contribution to adult sex ratio bias") +
  xlab("Sex bias in parameter")

Code
# Determine how much larger the contribution of each vital rates is compared to juvenile survival
# groundling vs nestling:
LTRE_coucal_ASR_prop <- 
  filter(LTRE_coucal_ASR) %>% 
  dplyr::select(-model) %>% 
  mutate(parameter2 = as.character(parameter)) %>% 
  dplyr::rename(parameter1 = parameter)

prop_contribution_table <- 
  LTRE_coucal_ASR_prop %>% 
  tidyr::expand(., parameter1 = parameter1, parameter2 = parameter1) %>%
  left_join(., dplyr::select(LTRE_coucal_ASR_prop, parameter1, contribution, sex, species), 
            by = c("parameter1")) %>% 
  dplyr::rename(contribution1 = contribution) %>% 
  left_join(., dplyr::select(LTRE_coucal_ASR_prop, parameter2, contribution, sex, species), 
            by = c("parameter2", "species", "sex")) %>% 
  dplyr::rename(contribution2 = contribution) %>% 
  mutate(prop_diff = abs(contribution1) / abs(contribution2),
         prop_diff2 = ifelse(contribution1 < 0 & contribution2 < 0, -(contribution1/contribution2),
                         ifelse(contribution1 < 0 & contribution2 > 0, contribution1/contribution2,
                                ifelse(contribution1 > 0 & contribution2 < 0, -(contribution1/contribution2), contribution1/contribution2))),
         parameter1 = factor(parameter1, 
                             levels = c("Hatching sex ratio",
                                        "Nestling survival",
                                        "Groundling survival",
                                        "Fledgling survival",
                                        "Adult survival",
                                        "Mating system")),
         parameter2 = factor(parameter2, 
                             levels = c("Hatching sex ratio",
                                        "Nestling survival",
                                        "Groundling survival",
                                        "Fledgling survival",
                                        "Adult survival",
                                        "Mating system")))

LTRE_heatmap <- 
  prop_contribution_table %>% 
  filter(parameter1 %in% c("Mating system",
                           "Adult survival",
                           "Fledgling survival",
                           "Groundling survival",
                           "Nestling survival",
                           "Hatching sex ratio") & 
                  parameter2 %in% c("Mating system",
                                    "Adult survival",
                                    "Fledgling survival",
                                    "Groundling survival",
                                    "Nestling survival",
                                    "Hatching sex ratio")) %>% 
  ggplot(., 
         aes(parameter1, parameter2, fill = prop_diff2)) + 
  geom_tile() +
  facet_grid(sex ~ species,
             labeller = labeller(.cols = species_names, .rows = analysis_names)) +
  theme(
    axis.text.x  = element_text(size = 8, angle = 45, hjust = 1),
    axis.text.y  = element_text(size = 8),
    axis.title.x  = element_text(size = 12),
    axis.title.y  = element_text(size = 12),
    legend.position = "bottom",
    legend.box = "horizontal",
    legend.text = element_text(size = 7),
    legend.title = element_text(size = 7),
    axis.ticks = element_blank(),
    panel.border = element_rect(color = "grey40", fill = NA),
    strip.background = element_blank(),
    strip.text = element_text(size = 12, colour = "black", face = "italic")) +
  scale_x_discrete(labels = c("Hatching sex ratio" = expression(italic("\u03C1")),
                              "Nestling survival" = expression(italic(S["n"])),
                              "Groundling survival" = expression(italic(S["g"])),
                              "Fledgling survival" = expression(italic(S["f"])),
                              "Adult survival" = expression(italic(phi["a"])),
                              "Mating system" = expression(italic("h")))) +
  scale_y_discrete(labels = c("Hatching sex ratio" = expression(italic("\u03C1")),
                              "Nestling survival" = expression(italic(S["n"])),
                              "Groundling survival" = expression(italic(S["g"])),
                              "Fledgling survival" = expression(italic(S["f"])),
                              "Adult survival" = expression(italic(phi["a"])),
                              "Mating system" = expression(italic("h")))) +
  scale_fill_gradient2(low = sex_pal2[1],
                       mid = "white",
                      high = sex_pal2[2], 
                      name = "Proportional contribution of A compared to B") +#,
  xlab("Parameter A") +
  ylab("Parameter B")

LTRE_heatmap

Supplementary Analysis 1: sex-specific juvenile movement

Here we assess individual heterogeneity in several movement traits based on radio telemetry data of tagged juveniles

Code
coucal_wp <-
  read.csv("data/raw/juvies_2013_2020_waypoints.csv",
           stringsAsFactors = FALSE, header = TRUE) %>%

  # clean up a few inconsistencies in the data
  mutate(timestamp = ifelse(ring_ID == "GN 76660" & date_dec == "20140502", "2014-05-02 10:09:00.000", timestamp),
         location.long = ifelse(location.long == 34.68357, 34.08357, location.long)) %>%
  filter(ring_ID != "GN 76726") %>%
  filter(ring_ID != "GN 76728") %>%
  filter(ring_ID != "GN 76732") %>%
  filter(ring_ID != "GN 76733") %>%
  filter(ring_ID != "GN 76777") %>% # both male and female??
  filter(movebank_event.id != "16897537877") %>%
  filter(movebank_event.id != "16897535798") %>%
  filter(movebank_event.id != "16897537758") %>%
  filter(movebank_event.id != "16897538205") %>%
  
  # specify the timezone of the time stamp
  mutate(timestamp = as.POSIXct(timestamp,
                                tz = "Africa/Nairobi"),

         # make the serial number a 4 digit string
         serial_number = as.character(str_pad(serial_number, 4, pad = "0"))) %>%
  
  # make a julian date time variable for plotting
  mutate(julian_date = as.numeric(format(timestamp, "%j"))) %>% 
  
  # remove observations older than 70 days of age
  filter(age_days < 71) %>%

  # consolidate columns of interest
  dplyr::select(species, year, nest.nr, ring_ID, sex, timestamp, julian_date, 
                age_days, location.long, location.lat, species, site, location, 
                comment)

# specify coordinate columns
coordinates(coucal_wp) <- c("location.long","location.lat")

# define spatial projection as WGS84
proj4string(coucal_wp) <- CRS("+init=epsg:4326")

Explore the spatial extent of the data

Code
# specify mapview options for viewing
mapviewOptions(basemaps = c("Esri.WorldImagery"),
               layers.control.pos = "topright",
               legend.pos = "bottomleft")

# check out the spatial data
mapview(coucal_wp[coucal_wp$ring_ID %in% "GN 76751",], zcol = "julian_date",
        layer.name = "WBC radio tagging",
        layers.control.pos = "topright")

Set-up the spatial data and create a trajectory format for analysis

Code
# make a spatial trajectory object of each animal's encounter history
tr_coucal_wp <- 
  as.ltraj(xy = sp::coordinates(coucal_wp),
           date = coucal_wp$timestamp,
           id = coucal_wp$ring_ID)

# convert the trajectory object back to a dataframe, rename columns, and
# calculate the cumulative distance traveled over the observation
tr_coucal_wp <- 
  ld(tr_coucal_wp) %>% 
  dplyr::rename(distance = dist,
                dispersion = R2n,
                cardinal_angle = abs.angle,
                relative_angle = rel.angle,
                ring_ID = id) %>% 
  group_by(ring_ID) %>% 
  mutate(cum_distance = cumsum(distance))

# extract the temporal summary info for each bird (min age, min date, max age) 
# and merge to the trajectory dataframe
coucal_wp_clean <- 
  as.data.frame(coucal_wp) %>% 
  group_by(species, ring_ID, sex, nest.nr) %>% 
  dplyr::summarise(min_age = min(age_days),
            min_date = min(timestamp),
            max_age = max(age_days)) %>% 
  arrange(desc(max_age)) %>% 
  left_join(., tr_coucal_wp, by = "ring_ID") %>% 
  
  # calcuate the differences in days between the observation and minimum date
  # to get time since first observation, convert from seconds to days, then add
  # to the minimum age value to extract the age at a given observation
  mutate(date_diff = date - min_date,
         year = year(date)) %>% 
  mutate(date_diff = round(((as.numeric(date_diff)/60)/60)/24)) %>% 
  mutate(age = min_age + date_diff) %>% 
  ungroup() %>% 
  
  # remove all NA rows
  filter(!is.na(cum_distance)) %>% 
  
  # calculate the temporal lag between observations
  group_by(ring_ID) %>% 
  mutate(obs_diff = date_diff - lag(date_diff),
         origin_x = x[which.min(date)],
         origin_y = y[which.min(date)])

# for each bird calculate the daily distance from origin
coucal_wp_clean$distance_origin <- 
  sapply(1:nrow(coucal_wp_clean), function(i)
    spDistsN1(as.matrix(coucal_wp_clean[i, c("x", "y")]), 
              as.matrix(coucal_wp_clean[i, c("origin_x", "origin_y")]), 
              longlat = TRUE))

# transform the distance metrics to meters
coucal_wp_clean <- 
  coucal_wp_clean %>% 
  mutate(distance = distance * 100000,
         distance_origin = distance_origin * 1000,
         cum_distance = cum_distance * 100000)

# check the distribution of distances
hist(log(coucal_wp_clean$distance))

Code
range(coucal_wp_clean$distance)
boxplot.stats(coucal_wp_clean$distance)
Code
# extract individuals that mainly had observations spaced 2 days apart
two_day_rings <- 
  coucal_wp_clean %>% 
  dplyr::select(species, ring_ID, sex, date_diff, obs_diff, age) %>% 
  group_by(ring_ID) %>% 
  dplyr::summarise(mean_obs_diff = mean(obs_diff, na.rm = TRUE),
            median_obs_diff = median(obs_diff, na.rm = TRUE)) %>% 
  arrange(desc(median_obs_diff)) %>% 
  filter(median_obs_diff <= 2)

# subset to only include individuals with mainly observations spaced 2 days apart
coucal_wp_clean <- 
  coucal_wp_clean %>% 
  filter(ring_ID %in% two_day_rings$ring_ID) %>% 
  
  # remove questionable observation 18 km 
  filter(distance < 18000)

Model 1: daily distance moved

Does daily distance moved increase with age since leaving the nest and does this relationship differ between males and females?

Code
#### Modeling sex-specific movements ####
# M1 Black Coucals: strong effect of age, no interaction with sex
mod_BC_dist <- 
  lmer(log(distance) ~ poly(age, 2) * sex + (1|ring_ID),
         data = filter(coucal_wp_clean, species == "BC"))

# effect-size plots
coefplot::coefplot(mod_BC_dist)

Code
plot(allEffects(mod_BC_dist))

Code
# extract fitted values
mod_BC_dist_fits <- 
  as.data.frame(effect(term = "poly(age, 2) * sex", mod = mod_BC_dist, 
                       xlevels = list(age = seq(min(coucal_wp_clean[coucal_wp_clean$species == "BC", "age"]),
                                                max(coucal_wp_clean[coucal_wp_clean$species == "BC", "age"]), 1),
                                      sex = c("Female", "Male"))))
mod_BC_dist_fits$species <- "BC"

# M1 White-browed Coucals: strong effect of age, no interaction with sex
mod_WBC_dist <- 
  lmer(log(distance) ~ poly(age, 2) * sex + 
         (1|ring_ID), 
       data = filter(coucal_wp_clean, species == "WBC" & distance > 0))

# effect-size plots
coefplot::coefplot(mod_WBC_dist)

Code
plot(allEffects(mod_WBC_dist))

Code
# extract fitted values
mod_WBC_dist_fits <- 
  as.data.frame(effect(term = "poly(age, 2) * sex", mod = mod_WBC_dist, 
                       xlevels = list(age = seq(min(coucal_wp_clean[coucal_wp_clean$species == "WBC", "age"]),
                                                max(coucal_wp_clean[coucal_wp_clean$species == "WBC", "age"]), 1),
                                      sex = c("Female", "Male"))))
mod_WBC_dist_fits$species <- "WBC"

mod_dist_fits <- 
  bind_rows(mod_BC_dist_fits, mod_WBC_dist_fits)
Code
# M1 plot
ggplot() +
  luke_theme +
  theme(legend.background = element_rect(fill="transparent")) +
  geom_line(data = coucal_wp_clean, 
            aes(x = age, y = log(distance), group = ring_ID, color = sex),
            alpha = 0.2) +
  geom_line(data = mod_dist_fits, aes(x = age, y = fit, color = sex),
            lwd = 0.5) +
  geom_ribbon(data = mod_dist_fits, aes(x = age, ymax = upper, ymin = lower, fill = sex),
              lwd = 1, alpha = 0.25) +
  facet_grid(species ~ ., 
             labeller = labeller(species = species.labs)) +
  ylab(expression(paste("Distance moved between observations (log10)" %+-%  "95% CI", sep = ""))) +
  xlab("Age since leaving nest") +
  scale_color_manual(values = sex_pal2,
                     name = "Sex",
                     breaks = c("Female", "Male"),
                     labels = c("Female", "Male")) +
  scale_fill_manual(values = sex_pal2,
                     name = "Sex",
                     breaks = c("Female", "Male"),
                     labels = c("Female", "Male"))

Code
# # first set up plotting parameters
# # define the plotting theme to be used in subsequent ggplots
# luke_theme <- 
#   theme_bw() +
#   theme(
#     text = element_text(family = "Verdana"),
#     legend.title = element_text(size = 16),
#     legend.text = element_text(size = 12),
#     axis.title.x = element_text(size = 12),
#     axis.text.x  = element_text(size = 12), 
#     axis.title.y = element_text(size = 16),
#     axis.text.y = element_text(size = 12),
#     strip.text = element_text(size = 12),
#     panel.grid.major = element_blank(),
#     panel.grid.minor = element_blank(),
#     axis.ticks = element_line(size = 0.5, colour = "grey40"),
#     axis.ticks.length = unit(0.2, "cm"),
#     panel.border = element_rect(linetype = "solid", colour = "grey"),
#     legend.position = c(0.1, 0.9)
#   )
# 
# # set plotting color palettes
# sex_pal2 <- RColorBrewer::brewer.pal(8, "Dark2")[c(2,1)]
# plot_palette_brood <- RColorBrewer::brewer.pal(7, "Blues")[c(3:7)]
# 
# # specify the facet labels for each species

Model 2: cumulative distance moved

Does cumulative distance moved increase with age since leaving the nest and does this relationship differ between males and females?

Code
# M2 Black Coucals: strong effect of age, no interaction with sex
mod_BC_cum_dist <- 
  lmer(log(cum_distance) ~ poly(age, 2) * sex + (1|ring_ID),
       data = filter(coucal_wp_clean, species == "BC"))

# effect-size plots
coefplot::coefplot(mod_BC_cum_dist)

Code
plot(allEffects(mod_BC_cum_dist))

Code
# extract fitted values
mod_BC_cum_dist_fits <- 
  as.data.frame(effect(term = "poly(age, 2) * sex", mod = mod_BC_cum_dist, 
                       xlevels = list(age = seq(min(coucal_wp_clean[coucal_wp_clean$species == "BC", "age"]),
                                                max(coucal_wp_clean[coucal_wp_clean$species == "BC", "age"]), 1),
                                      sex = c("Female", "Male"))))
mod_BC_cum_dist_fits$species <- "BC"

# M2 White-browed Coucals: strong effect of age, no interaction with sex
mod_WBC_cum_dist <- 
  lmer(log(cum_distance) ~ poly(age, 2) * sex + 
         (1|ring_ID), 
       data = filter(coucal_wp_clean, species == "WBC" & distance > 0))

# effect-size plots
coefplot::coefplot(mod_WBC_cum_dist)

Code
plot(allEffects(mod_WBC_cum_dist))

Code
# extract fitted values
mod_WBC_cum_dist_fits <- 
  as.data.frame(effect(term = "poly(age, 2) * sex", mod = mod_WBC_cum_dist, 
                       xlevels = list(age = seq(min(coucal_wp_clean[coucal_wp_clean$species == "WBC", "age"]),
                                                max(coucal_wp_clean[coucal_wp_clean$species == "WBC", "age"]), 1),
                                      sex = c("Female", "Male"))))
mod_WBC_cum_dist_fits$species <- "WBC"

mod_cum_dist_fits <- 
  bind_rows(mod_BC_cum_dist_fits, mod_WBC_cum_dist_fits)
Code
# M2 plot
ggplot() +
  luke_theme +
  theme(legend.background = element_rect(fill="transparent")) +
  geom_line(data = coucal_wp_clean, 
            aes(x = age, y = log(cum_distance), group = ring_ID, color = sex),
            alpha = 0.2) +
  geom_line(data = mod_cum_dist_fits, aes(x = age, y = fit, color = sex),
            lwd = 0.5) +
  geom_ribbon(data = mod_cum_dist_fits, aes(x = age, ymax = upper, ymin = lower, fill = sex),
              lwd = 1, alpha = 0.25) +
  facet_grid(species ~ ., 
             labeller = labeller(species = species.labs)) +
  ylab(expression(paste("Cumulative distance moved (log10)" %+-%  "95% CI", sep = ""))) +
  xlab("Age since leaving nest") +
  scale_color_manual(values = sex_pal2,
                     name = "Sex",
                     breaks = c("Female", "Male"),
                     labels = c("Female", "Male")) +
  scale_fill_manual(values = sex_pal2,
                    name = "Sex",
                    breaks = c("Female", "Male"),
                    labels = c("Female", "Male"))

Model 3: distance from origin

Do individuals move further from the nest as they age and does this relationship differ between males and females?

Code
# M3 Black Coucals: strong effect of age, no interaction with sex
mod_BC_disp <- 
  lmer(log(distance_origin) ~ poly(age, 2) * sex + (1|ring_ID),
       data = filter(coucal_wp_clean, species == "BC" & distance_origin > 0))

# effect-size plots
coefplot::coefplot(mod_BC_disp)

Code
plot(allEffects(mod_BC_disp))

Code
# extract fitted values
mod_BC_disp_fits <- 
  as.data.frame(effect(term = "poly(age, 2) * sex", mod = mod_BC_disp, 
                       xlevels = list(age = seq(min(coucal_wp_clean[coucal_wp_clean$species == "BC", "age"]),
                                                max(coucal_wp_clean[coucal_wp_clean$species == "BC", "age"]), 1),
                                      sex = c("Female", "Male"))))
mod_BC_disp_fits$species <- "BC"

# M3 White-browed Coucals: strong effect of age, no interaction with sex
mod_WBC_disp <- 
  lmer(log(distance_origin) ~ poly(age, 2) * sex + 
         (1|ring_ID), 
       data = filter(coucal_wp_clean, species == "WBC" & distance_origin > 0))

# effect-size plots
coefplot::coefplot(mod_WBC_disp)

Code
plot(allEffects(mod_WBC_disp))

Code
# extract fitted values
mod_WBC_disp_fits <- 
  as.data.frame(effect(term = "poly(age, 2) * sex", mod = mod_WBC_disp, 
                       xlevels = list(age = seq(min(coucal_wp_clean[coucal_wp_clean$species == "WBC", "age"]),
                                                max(coucal_wp_clean[coucal_wp_clean$species == "WBC", "age"]), 1),
                                      sex = c("Female", "Male"))))
mod_WBC_disp_fits$species <- "WBC"

mod_disp_fits <- 
  bind_rows(mod_BC_disp_fits, mod_WBC_disp_fits)
Code
# M3 plot
ggplot() +
  luke_theme +
  theme(legend.background = element_rect(fill="transparent")) +
  geom_line(data = filter(coucal_wp_clean, distance_origin > 0), 
            aes(x = age, y = log(distance_origin), group = ring_ID, color = sex),
            alpha = 0.2) +
  geom_line(data = mod_disp_fits, aes(x = age, y = fit, color = sex),
            lwd = 0.5) +
  geom_ribbon(data = mod_disp_fits, aes(x = age, ymax = upper, ymin = lower, fill = sex),
              lwd = 1, alpha = 0.25) +
  facet_grid(species ~ ., 
             labeller = labeller(species = species.labs)) +
  ylab(expression(paste("Daily distance from nest (log10)" %+-%  "95% CI", sep = ""))) +
  xlab("Age since leaving nest") +
  scale_color_manual(values = sex_pal2,
                     name = "Sex",
                     breaks = c("Female", "Male"),
                     labels = c("Female", "Male")) +
  scale_fill_manual(values = sex_pal2,
                    name = "Sex",
                    breaks = c("Female", "Male"),
                    labels = c("Female", "Male"))

Supplementary Analysis 2: weekly sex-specific adult survival

Here we run a weekly mark-recapture survival analysis on the weekly encounters (and dead recoveries) of adults throughout the season. The first step is to wrangle the capture history into the format needed for a Burnham-joint-live-dead model: each encounter event requires two characters (00 = not encountered, 01 = encountered alive, 11 = encountered dead)

Code
status_dat_ad <- 
  # read raw data
  read.csv("data/raw/Coucal_adults_survival_2001-2019_20200129.csv", 
           header = TRUE, stringsAsFactors = FALSE, na.strings = c("", " ", "NA")) %>%
  
  # rename ring_ID column
  dplyr::rename(ring_ID = Alu,
                ageC2 = days_from_ringing_until_last_observation) %>% 
  
  # specify date strings properly
  mutate(month_ringed = str_sub(date_dec_initially_ringed, start = 5, end = 6),
         day_ringed = str_sub(date_dec_initially_ringed, start = 7, end = 8),
         year_ringed = str_sub(date_dec_initially_ringed, start = 1, end = 4),
         month_last_ob = str_sub(date_dec_last_observed, start = 5, end = 6),
         day_last_ob = str_sub(date_dec_last_observed, start = 7, end = 8),
         year_last_ob = str_sub(date_dec_initially_ringed, start = 1, end = 4)) %>% 
  mutate(date_ringed = as.Date(paste(year_ringed, month_ringed, day_ringed, sep = "-"), 
                               format = "%Y-%m-%d"),
         date_last_ob = as.Date(paste(year_last_ob, month_last_ob, day_last_ob, sep = "-"), 
                                format = "%Y-%m-%d")) %>% 
  
  # select variables of interest
  dplyr::select(species, ring_ID, sex, year, ageC2, ageC, statusC, date_ringed, date_last_ob) %>% 
  
  # remove all white space from data
  mutate(across(everything(), ~str_trim(.x))) %>% 
  mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>% 
  
  # specify empty data as NA
  mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>% 
  
  # exclude all individuals that died in the nest
  # filter(Fledged_status == "yes") %>% 
  
  # classify columns
  mutate(sex = as.factor(sex),
         ageC2 = as.numeric(ageC2)) %>% 
  
  # remove rows with missing sex, age, and status info
  filter(!is.na(sex) & !is.na(ageC2) & !is.na(statusC)) %>% 
  
  # make a unique id for each individual
  mutate(# create the age of entry into the data (all at age 15)
    entry = 0,
    
    # specify the age of death or censoring
    exit = ageC2,
    
    # make the event numeric and specify if 
    # the individual died (1) or was censored (0)
    event = as.numeric(statusC),
    
    sex = ifelse(str_detect(sex, pattern = "[Ff]emale"), "F", 
                 ifelse(str_detect(sex, pattern = "[Mm]ale"), "M", sex))) %>% 
  
  # consolidate to variables of interest
  dplyr::select(species, ring_ID, year, sex, entry, exit, event, date_ringed, date_last_ob)

status_dat_ad_week <- 
  status_dat_ad %>% 
  dplyr::select(species, ring_ID, sex, date_ringed, date_last_ob) %>% 
  pivot_longer(-c(species, ring_ID, sex), names_to = "status", values_to = "date") %>%
  mutate(year_numeric = as.numeric(lubridate::year(date)),
         week_numeric = as.numeric(strftime(date, format = "%V"))) %>% 
  dplyr::select(species, ring_ID, year_numeric, sex, week_numeric)

# import raw csv data into R
BC_detect_dat_A <- 
  read_xlsx("data/raw/All_coucal_waypoints_2001_2019_20200202.xlsx", na = "NA", col_types = "text") %>% 
  dplyr::select(species, ring_ID, sex, year, site, age_status, date_dec) %>% 
  mutate(month = str_sub(date_dec, start = 5, end = 6),
         day = str_sub(date_dec, start = 7, end = 8)) %>% 
  mutate(date = as.Date(paste(year, month, day, sep = "-"), format = "%Y-%m-%d")) %>% 
  mutate(across(c("sex", "site", "age_status"), tolower)) %>%
  mutate(age_status = ifelse(age_status == "adult", "A", ifelse(age_status == "juvenile", "J", "XXX")),
         ring_ID = str_replace_all(string = ring_ID, fixed(" "), "")) %>% 
  mutate(across(c("species", "ring_ID", "sex", "site", "age_status"), as.factor)) %>% 
  mutate(sex = ifelse(str_detect(sex, pattern = "[Ff]emale"), "F",
                      ifelse(str_detect(sex, pattern = "[Mm]ale"), "M", sex)),
         month_year = format(date, "%Y-%m"),
         month_numeric = as.numeric(month),
         year_numeric = as.numeric(year),
         week_numeric = as.numeric(strftime(date, format = "%V"))) %>% 
  filter(species == "BC") %>% 
  filter(age_status == "A") %>% 
  dplyr::select(species, ring_ID, year_numeric, sex, week_numeric) %>%
  bind_rows(., filter(status_dat_ad_week, species == "BC")) %>% 
  filter(!is.na(week_numeric)) %>% 
  distinct() %>% 
  mutate(week_std = round(scale_by(week_numeric ~ year_numeric, ., scale = 0), digits = 0)) %>% 
  mutate(week_std = week_std + abs(min(week_std, na.rm = TRUE))) %>% 
  dplyr::rename(year = year_numeric)

WBC_detect_dat_A <- 
  read_xlsx("data/raw/All_coucal_waypoints_2001_2019_20200202.xlsx", na = "NA", col_types = "text") %>% 
  dplyr::select(species, ring_ID, sex, year, site, age_status, date_dec) %>% 
  mutate(month = str_sub(date_dec, start = 5, end = 6),
         day = str_sub(date_dec, start = 7, end = 8)) %>% 
  mutate(date = as.Date(paste(year, month, day, sep = "-"), format = "%Y-%m-%d")) %>% 
  mutate(across(c("sex", "site", "age_status"), tolower)) %>%
  mutate(age_status = ifelse(age_status == "adult", "A", ifelse(age_status == "juvenile", "J", "XXX")),
         ring_ID = str_replace_all(string = ring_ID, fixed(" "), "")) %>% 
  mutate(across(c("species", "ring_ID", "sex", "site", "age_status"), as.factor)) %>% 
  mutate(sex = ifelse(str_detect(sex, pattern = "[Ff]emale"), "F",
                      ifelse(str_detect(sex, pattern = "[Mm]ale"), "M", sex)),
         month_year = format(date, "%Y-%m"),
         month_numeric = as.numeric(month),
         year_numeric = as.numeric(year),
         week_numeric = as.numeric(strftime(date, format = "%V"))) %>% 
  filter(species == "WBC") %>% 
  filter(age_status == "A") %>% 
  dplyr::select(species, ring_ID, year_numeric, sex, week_numeric) %>%
  bind_rows(., filter(status_dat_ad_week, species == "WBC")) %>% 
  filter(!is.na(week_numeric)) %>% 
  distinct() %>% 
  mutate(week_std = round(scale_by(week_numeric ~ year_numeric, ., scale = 0), digits = 0)) %>% 
  mutate(week_std = week_std + abs(min(week_std, na.rm = TRUE))) %>% 
  dplyr::rename(year = year_numeric)

# use the BaSTA function "CensusToCaptHist()" function to convert long format
# encounter histories of each individual, to wide format with 1's and 0's for 
# each week of encounter for each year
BC_weekly_A_ch <- 
  CensusToCaptHist(ID = BC_detect_dat_A$ring_ID,
                   d = BC_detect_dat_A$week_std,
                   dformat = "%W", timeInt = "%W") %>% 
  dplyr::rename(ring_ID = ID) %>%
  mutate(ID = rownames(.)) %>% 
  left_join(., dplyr::select(BC_detect_dat_A, ring_ID, sex, year), by = "ring_ID") %>% 
  distinct()

WBC_weekly_A_ch <- 
  CensusToCaptHist(ID = WBC_detect_dat_A$ring_ID,
                   d = WBC_detect_dat_A$week_std,
                   dformat = "%W", timeInt = "%W") %>% 
  dplyr::rename(ring_ID = ID) %>%
  mutate(ID = rownames(.)) %>% 
  left_join(., dplyr::select(WBC_detect_dat_A, ring_ID, sex, year), by = "ring_ID") %>% 
  distinct()

# create a two-character string for each encounter and clean the output
BC_weekly_A_ch_Burnham <- 
  sapply(BC_weekly_A_ch[2:26], function(x) paste(x, "0", sep = "")) %>% 
  cbind(BC_weekly_A_ch[c(1,45:length(BC_weekly_A_ch))]) %>% 
  mutate(across(everything(), ~str_replace(string = .x, pattern = "NA0", "00"))) %>% 
  mutate(across(everything(), ~str_trim(.x))) %>% 
  mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>% 
  mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>% 
  distinct() %>% 
  mutate(sex = as.factor(sex))

WBC_weekly_A_ch_Burnham <- 
  sapply(WBC_weekly_A_ch[2:26], function(x) paste(x, "0", sep = "")) %>% 
  cbind(WBC_weekly_A_ch[c(1,54:length(WBC_weekly_A_ch))]) %>% 
  mutate(across(everything(), ~str_replace(string = .x, pattern = "NA0", "00"))) %>% 
  mutate(across(everything(), ~str_trim(.x))) %>% 
  mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>% 
  mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>% 
  distinct() %>% 
  mutate(sex = as.factor(sex))

# Import status data
status_dat_all_ad <- 
  # read raw data
  read.csv("data/raw/Coucal_adults_survival_2001-2019_20200129.csv", 
           header = TRUE, stringsAsFactors = FALSE, na.strings = c("", " ", "NA")) %>%
  
  # rename ring_ID column
  dplyr::rename(ring_ID = Alu,
                ageC2 = days_from_ringing_until_last_observation) %>% 
  
  # specify date strings properly
  mutate(month_ringed = str_sub(date_dec_initially_ringed, start = 5, end = 6),
         day_ringed = str_sub(date_dec_initially_ringed, start = 7, end = 8),
         year_ringed = str_sub(date_dec_initially_ringed, start = 1, end = 4),
         month_last_ob = str_sub(date_dec_last_observed, start = 5, end = 6),
         day_last_ob = str_sub(date_dec_last_observed, start = 7, end = 8),
         year_last_ob = str_sub(date_dec_initially_ringed, start = 1, end = 4)) %>% 
  mutate(date_ringed = as.Date(paste(year_ringed, month_ringed, day_ringed, sep = "-"), 
                               format = "%Y-%m-%d"),
         date_last_ob = as.Date(paste(year_last_ob, month_last_ob, day_last_ob, sep = "-"), 
                                format = "%Y-%m-%d")) %>% 
  
  # select variables of interest
  dplyr::select(species, ring_ID, sex, year, ageC2, ageC, statusC, date_ringed, date_last_ob) %>% 
  
  # remove all white space from data
  mutate(across(everything(), ~str_trim(.x))) %>% 
  mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>% 
  
  # specify empty data as NA
  mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>% 
  
  # classify columns
  mutate(sex = as.factor(sex),
         ageC2 = as.numeric(ageC2)) %>% 
  
  # remove rows with missing sex, age, and status info
  filter(!is.na(sex) & !is.na(ageC2) & !is.na(statusC)) %>% 
  
  # make a unique id for each individual
  mutate(# create the age of entry into the data (all at age 15)
    entry = 0,
    
    # specify the age of death or censoring
    exit = ageC2,
    
    # make the event numeric and specify if 
    # the individual died (1) or was censored (0)
    event = as.numeric(statusC),
    
    sex = ifelse(str_detect(sex, pattern = "[Ff]emale"), "F", 
                 ifelse(str_detect(sex, pattern = "[Mm]ale"), "M", sex))) %>% 
  
  # consolidate to variables of interest
  dplyr::select(species, ring_ID, year, sex, entry, exit, event, date_ringed, date_last_ob) %>% 
  dplyr::select(species, ring_ID, year, sex, event)

# join detection history with status data
BC_detect_status_join <-
  filter(status_dat_all_ad, species == "BC") %>%
  left_join(BC_weekly_A_ch_Burnham, ., by = c("ring_ID", "year", "sex")) %>% 
  dplyr::select(-species) %>% 
  mutate(event = ifelse(is.na(event), 0, event))

WBC_detect_status_join <-
  filter(status_dat_all_ad, species == "WBC") %>%
  left_join(WBC_weekly_A_ch_Burnham, ., by = c("ring_ID", "year", "sex")) %>% 
  dplyr::select(-species) %>% 
  mutate(event = ifelse(is.na(event), 0, event))

# determine the last and first detection for each individual 
max_age_index_BC <- 
  apply(BC_detect_status_join[, c(1:25)], 1, function(x) which(x == "10")) %>% 
  lapply(., function(x) x[which.max(x)]) %>% 
  unlist(.)

min_age_index_BC <- 
  apply(BC_detect_status_join[, c(1:25)], 1, function(x) which(x == "10")) %>% 
  lapply(., function(x) x[which.min(x)]) %>% 
  unlist(.)

max_age_index_WBC <- 
  apply(WBC_detect_status_join[, c(1:25)], 1, function(x) which(x == "10")) %>% 
  lapply(., function(x) x[which.max(x)]) %>% 
  unlist(.)

min_age_index_WBC <- 
  apply(WBC_detect_status_join[, c(1:25)], 1, function(x) which(x == "10")) %>% 
  lapply(., function(x) x[which.min(x)]) %>% 
  unlist(.)

# put "11" at the last detection if the status was a "1" (i.e., dead)
for(i in 1:nrow(BC_detect_status_join)){
  if(BC_detect_status_join$event[i] == "1"){
    BC_detect_status_join[i, max_age_index_BC[i]] <- "11"
  }
}

for(i in 1:nrow(WBC_detect_status_join)){
  if(WBC_detect_status_join$event[i] == "1"){
    WBC_detect_status_join[i, max_age_index_WBC[i]] <- "11"
  }
}

# bind the min and max age info to the capture history data
detect_status_BC <-
  bind_cols(BC_detect_status_join, max_age_index_BC) %>% 
  bind_cols(., min_age_index_BC) %>% 
  dplyr::rename(std_week_last_obs = '...31',
                std_week_first_obs = '...32')

detect_status_WBC <-
  bind_cols(WBC_detect_status_join, max_age_index_WBC) %>% 
  bind_cols(., min_age_index_WBC) %>% 
  dplyr::rename(std_week_last_obs = '...31',
                std_week_first_obs = '...32')

# consolidate capture histories for each sex
Black_Coucal_adult_weekly_Burnham_ch <-
  apply(detect_status_BC[, c(1:25)], 1, paste, collapse = "") %>% 
  as.data.frame() %>% 
  dplyr::rename(ch = '.') %>% 
  mutate(ch = as.character(ch)) %>% 
  dplyr::bind_cols(., detect_status_BC[, c(26:length(detect_status_BC))]) %>% 
  filter(str_detect(ch, "11", negate = TRUE) | str_count(ch, "1") > 2) %>%
  dplyr::select(ch, ring_ID, sex, year) %>%
  drop_na() %>% 
  mutate(year = ifelse(year == "1006", "2006", year)) %>%
  mutate_at(vars(ring_ID, sex, year), as.factor) 

White_browed_Coucal_adult_weekly_Burnham_ch <-
  apply(detect_status_WBC[, c(1:25)], 1, paste, collapse = "") %>% 
  as.data.frame() %>% 
  dplyr::rename(ch = '.') %>% 
  mutate(ch = as.character(ch)) %>% 
  dplyr::bind_cols(., detect_status_WBC[, c(26:length(detect_status_WBC))]) %>% 
  filter(str_detect(ch, "11", negate = TRUE) | str_count(ch, "1") > 2) %>%
  dplyr::select(ch, ring_ID, sex, year) %>%
  drop_na() %>%
  mutate_at(vars(ring_ID, sex, year), as.factor)

BC_adult_ch <-
  Black_Coucal_adult_weekly_Burnham_ch %>%
  mutate(sex = as.factor(sex),
         year = as.factor(year)) %>%
  dplyr::select(ch, sex, year)

WBC_adult_ch <-
  White_browed_Coucal_adult_weekly_Burnham_ch %>%
  mutate(sex = as.factor(sex),
         year = as.factor(year)) %>%
  dplyr::select(ch, sex, year)

Run the survival analysis in RMark, extract the model predictions, and plot the results.

Code
# process the adult data as a "Burnham" analysis
BC_coucal_adult.proc <- RMark::process.data(data = BC_adult_ch,
                                                model = "Burnham",
                                                groups = c("sex", "year"))

WBC_coucal_adult.proc <- RMark::process.data(data = WBC_adult_ch,
                                             model = "Burnham",
                                             groups = c("sex", "year"))

# make design matrix
BC_coucal_adult.ddl <- 
  RMark::make.design.data(BC_coucal_adult.proc)

WBC_coucal_adult.ddl <- 
  RMark::make.design.data(WBC_coucal_adult.proc)

burnham_S_analysis_run_pFr_dot <- 
  function(proc_data, design_data){
    
    # sex-specific model
    S.sex = list(formula = ~sex)
    
    # sex- and week-specific model
    S.sex_Time = list(formula = ~sex + Time)
    
    # sex- and week-specific model (quadratic)
    S.sex_Quad = list(formula = ~sex + (I(Time) + I(Time^2)))
    
    # sex- and week-specific model (interaction)
    S.sex_x_Time = list(formula = ~sex * Time)
    
    # sex- and week-specific model (quadratic with interaction)
    S.sex_x_Quad = list(formula = ~sex * (I(Time) + I(Time^2)))
    
    # p (encounter probability):
    # constant model
    p.dot = list(formula = ~1)
    
    # F (site fidelity probability):
    # constant model
    F.dot = list(formula = ~1)
    
    # r (recovery probability)
    # constant model
    r.dot = list(formula = ~1)
    
    # create model list for all a priori models above
    cml <- RMark::create.model.list("Burnham")
    
    # run the models in program MARK
    model.list <-  RMark::mark.wrapper(cml, data = proc_data, 
                                       ddl = design_data, delete = TRUE, 
                                       wrap = FALSE, threads = 1, brief = TRUE,
                                       silent = TRUE, output = FALSE)
    
    # output the model list and sotre the results
    return(model.list)
  }

adult_S_analysis_out_BC <-
  burnham_S_analysis_run_pFr_dot(proc_data = BC_coucal_adult.proc,
                                 design_data = BC_coucal_adult.ddl)

adult_S_analysis_out_WBC <-
  burnham_S_analysis_run_pFr_dot(proc_data = WBC_coucal_adult.proc,
                                 design_data = WBC_coucal_adult.ddl)

BC_sex_x_quad_reals <- 
  adult_S_analysis_out_BC$S.sex_x_Quad.p.dot.r.dot.F.dot$results$real %>% 
  bind_cols(data.frame(str_split_fixed(rownames(.), " ", 
                                       n = 5)), .) %>% 
  mutate(age = as.integer(unlist(str_extract_all(X4,"[0-9]+"))),
         sex = as.factor(ifelse(unlist(str_extract_all(X2,"[FM]")) == "F", 
                                "Female", "Male"))) %>% 
  filter(X1 == "S") %>% 
  dplyr::select(sex, age, estimate, se, lcl, ucl) %>% 
  mutate(species = "Black Coucal")

WBC_sex_x_quad_reals <- 
  adult_S_analysis_out_WBC$S.sex_x_Quad.p.dot.r.dot.F.dot$results$real %>% 
  bind_cols(data.frame(str_split_fixed(rownames(.), " ", 
                                       n = 5)), .) %>% 
  mutate(age = as.integer(unlist(str_extract_all(X4,"[0-9]+"))),
         sex = as.factor(ifelse(unlist(str_extract_all(X2,"[FM]")) == "F", 
                                "Female", "Male"))) %>% 
  filter(X1 == "S") %>% 
  dplyr::select(sex, age, estimate, se, lcl, ucl) %>% 
  mutate(species = "White-browed Coucal")

BC_sex_reals <- 
  adult_S_analysis_out_BC$S.sex.p.dot.r.dot.F.dot$results$real %>% 
  bind_cols(data.frame(str_split_fixed(rownames(.), " ", 
                                       n = 5)), .) %>% 
  mutate(age = as.integer(unlist(str_extract_all(X4,"[0-9]+"))),
         sex = as.factor(ifelse(unlist(str_extract_all(X2,"[FM]")) == "F", 
                                "Female", "Male"))) %>% 
  filter(X1 == "S") %>% 
  dplyr::select(sex, age, estimate, se, lcl, ucl) %>% 
  mutate(species = "Black Coucal")

WBC_sex_reals <- 
  adult_S_analysis_out_WBC$S.sex.p.dot.r.dot.F.dot$results$real %>% 
  bind_cols(data.frame(str_split_fixed(rownames(.), " ", 
                                       n = 5)), .) %>% 
  mutate(age = as.integer(unlist(str_extract_all(X4,"[0-9]+"))),
         sex = as.factor(ifelse(unlist(str_extract_all(X2,"[FM]")) == "F", 
                                "Female", "Male"))) %>% 
  filter(X1 == "S") %>% 
  dplyr::select(sex, age, estimate, se, lcl, ucl) %>% 
  mutate(species = "White-browed Coucal")

time_varying_plot <- 
  ggplot(data = bind_rows(BC_sex_x_quad_reals, WBC_sex_x_quad_reals)) +
    geom_line(aes(x = age, y = estimate, group = sex, color = sex), lwd = 0.8) +
    geom_ribbon(aes(x = age, ymin = lcl, ymax = ucl, group = sex, fill = sex), alpha = 0.5) +
    scale_color_manual(values = sex_pal2,
                       labels = c("Female", "Male")) +
    scale_fill_manual(values = sex_pal2,
                      labels = c("Female", "Male")) +
  facet_grid(. ~ species) +
  theme_bw() +
  theme(text = element_text(family = "Verdana"),
        axis.title = element_text(size = 12),
        axis.text  = element_text(size = 10), 
        strip.text = element_text(size = 12),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks = element_line(size = 0.5, colour = "grey40"),
        axis.ticks.length = unit(0.2, "cm"),
        panel.border = element_rect(linetype = "solid", colour = "grey"),
        legend.position = c(0.8, 0.2),
        legend.title = element_blank()) +
  xlab("Week of breeding season (standardized across years)") +
  ylab("Adult weekly survival probability (± 95% CI)")

constant_plot <- 
  ggplot2::ggplot() +
    geom_errorbar(data = bind_rows(BC_sex_reals, WBC_sex_reals),
                  aes(x = sex, ymax = ucl, ymin = lcl, color = sex),
                  alpha = 1, width = 0.05, lwd = 0.5) +
    geom_point(data = bind_rows(BC_sex_reals, WBC_sex_reals),
               aes(x = sex, y = estimate, fill = sex),
               lwd = 3, shape = 21) +
    scale_color_manual(values = sex_pal2,
                       labels = c("Female", "Male")) +
    scale_fill_manual(values = sex_pal2,
                      labels = c("Female", "Male")) +
    facet_grid(. ~ species) +
    theme_bw() +
    theme(
      text = element_text(family = "Verdana"),
      axis.title = element_blank(),
      axis.text  = element_text(size = 10), 
      strip.text = element_text(size = 12),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      axis.ticks = element_line(size = 0.5, colour = "grey40"),
      axis.ticks.length = unit(0.2, "cm"),
      panel.border = element_rect(linetype = "solid", colour = "grey"),
      legend.position = "none") +
    ylab("Adult weekly survival probability (± 95% CI)") +
    scale_y_continuous(limits = c(0, 1))

combo_plot_S <- 
  time_varying_plot + constant_plot + plot_annotation(tag_levels = "A")
combo_plot_S

Supplementary Analysis 3: Plumage-based sex-specific age structure

Here we assess annual variation in the sex-specific age structure of the Black Coucal population. Older individuals (i.e., 3+ years) have rufous flight feathers (i.e., primary coverts, primaries, secondary coverts, secondaries), whereas these flight feathers are barred for younger individuals

Code
age_plumage_dist_data <-
  read_xlsx("data/raw/age_classes_2001-2020_2.xlsx", na = "NA", col_types = "text") %>%
  dplyr::rename_all(~str_replace_all(., "\\s+", "")) %>%
  dplyr::mutate(species = tolower(species)) %>%
  dplyr::filter(str_detect(string = species, pattern = "black")) %>% 
  dplyr::mutate(month = str_sub(`date(yymmdd)`, start = 5, end = 6),
         day = str_sub(`date(yymmdd)`, start = 7, end = 8)) %>% 
  dplyr::mutate(date = as.Date(paste(year, month, day, sep = "-"), format = "%Y-%m-%d")) %>% 
  # remove all white space from data
  dplyr::mutate(dplyr::across(.cols = dplyr::everything(), 
                str_remove_all, pattern = fixed(" "))) %>% 
  dplyr::mutate(species = "BC",
                sex = ifelse(sex == "female", "F", ifelse(sex == "male", "M", "XXX"))) %>% 
  dplyr::mutate(dplyr::across(.cols = dplyr::everything(), 
                              str_replace_all, pattern = fixed("onebarred"), replacement = fixed("barred"))) %>%
  dplyr::mutate(Secondarycoverts = ifelse(SN == "462", "rufous", Secondarycoverts)) %>% 
  dplyr::mutate(Secondarycoverts = ifelse(SN == "409", "barred", Secondarycoverts)) %>% 
  dplyr::mutate(age_LEP = ifelse(str_detect(primarycoverts, "barred") & 
                                   str_detect(primaries, "barred") & 
                                   str_detect(Secondarycoverts, "barred") & 
                                   str_detect(secondaries, "barred"), 1,
                                 
                                 ifelse(str_detect(primaries, "barred") | 
                                        str_detect(secondaries, "barred"), 2,
                                        
                                        ifelse((str_detect(primarycoverts, "barred") |
                                               str_detect(Secondarycoverts, "barred")) &
                                               (str_detect(primaries, "rufous") |
                                               str_detect(secondaries, "rufous")), 2,
                                               
                                               ifelse(str_detect(primarycoverts, "rufous") & 
                                                         str_detect(primaries, "rufous") & 
                                                         str_detect(Secondarycoverts, "rufous") & 
                                                         str_detect(secondaries, "rufous"), 3, 0))))) %>%
  dplyr::mutate(age = ifelse(age == "juvenile", 1, age)) %>%
  dplyr::mutate(CHECK = ifelse(age != age_LEP, "CHECK", "")) %>% 
  dplyr::mutate(age_LEP = age) %>% 
  dplyr::filter(!is.na(age_LEP)) %>% 
  dplyr::filter(year != "2020") %>% 
  dplyr::group_by(year, sex, age_LEP) %>%
  dplyr::summarise(n_ages = n_distinct(Alu)) %>% 
  dplyr::group_by(year, sex) %>% 
  dplyr::mutate(prop = n_ages/sum(n_ages)) %>%
  dplyr::mutate(age_LEP = as.character(age_LEP)) %>% 
  dplyr::mutate(age_LEP = ifelse(age_LEP == "3", "3+", age_LEP)) %>% 
  dplyr::mutate(age_LEP = fct_rev(as.factor(age_LEP)))

age_plumage_dist_data2 <- 
  age_plumage_dist_data %>% 
  dplyr::group_by(sex, age_LEP) %>% 
  dplyr::summarise(mean_n_ages = mean(n_ages)) %>% 
  dplyr::group_by(sex) %>% 
  dplyr::mutate(prop = mean_n_ages/sum(mean_n_ages),
                year = "Grand\naverage") %>% 
  dplyr::bind_rows(age_plumage_dist_data, .) %>% 
  dplyr::mutate(sex = ifelse(sex == "F", "Female", "Male"))

age_structure_plot <- 
  ggplot(age_plumage_dist_data2, 
       aes(fill = sex, alpha = age_LEP, 
           y = prop, x = sex)) + 
  geom_bar(position="stack", stat="identity") +
  facet_wrap(. ~ year) +
  theme_bw() +
  theme(text = element_text(family = "Verdana"),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.ticks.length = unit(0.1, "cm"),
    panel.border = element_blank(),
    panel.spacing = unit(0.3, "lines"),
    strip.background = element_blank(),
    legend.position = "top") +
  scale_fill_manual(values = sex_pal2, name = "Sex", guide = FALSE) +
  scale_alpha_manual(values = c(0.3, 0.65, 1), name = "Age (years)",
                     guide = guide_legend(title.position = "top", title.hjust = 0.5, label.position = "bottom")) +
  ylab("Proportion of breeding adults") +
  xlab("Sex")
age_structure_plot

Code
age_structure_table <- 
  age_plumage_dist_data2 %>% 
  dplyr::filter(str_detect(year, "Grand")) %>% 
  dplyr::ungroup() %>% 
  dplyr::select(sex, age_LEP, mean_n_ages, prop) %>% 
  dplyr::arrange(sex, desc(age_LEP)) %>% 
  gt(rowname_col = "row",
     groupname_col = "sex") %>%
  cols_label(age_LEP = html("<i>Age class</i>"), 
             mean_n_ages = "Annual\nmean no. ind.",
             prop = "Proportion") %>% 
  fmt_number(columns = vars(mean_n_ages, prop),
             decimals = 2,
             use_seps = FALSE) %>% 
  cols_align(align = "right",
             columns = vars(mean_n_ages)) %>%
  cols_align(align = "left",
             columns = vars(prop)) %>%
  tab_options(row_group.font.weight = "bold",
              row_group.background.color = brewer.pal(9,"Greys")[3],
              table.font.size = 12,
              data_row.padding = 3,
              row_group.padding = 4,
              summary_row.padding = 2,
              column_labels.font.size = 14,
              row_group.font.size = 12,
              table.width = pct(80))
age_structure_table
Age class Annual mean no. ind. Proportion
Female
1 2.38 0.21
2 3.67 0.32
3+ 5.36 0.47
Male
1 3.46 0.34
2 3.31 0.33
3+ 3.40 0.33

Housekeeping

Display R version and package information for time-dependent reproducibility

Code
sessionInfo()
R version 4.2.2 (2022-10-31)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Monterey 12.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  grid      stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] remotes_2.4.2.1      merDeriv_0.2-4       lavaan_0.6-16       
 [4] sandwich_3.0-2       nonnest2_0.5-6       mapview_2.11.0      
 [7] coefplot_1.2.8       reshape_0.8.9        adehabitatLT_0.3.27 
[10] CircStats_0.2-6      boot_1.3-28.1        adehabitatMA_0.3.16 
[13] ade4_1.7-22          sp_2.1-0             reshape2_1.4.4      
[16] magrittr_2.0.3       epitools_0.5-10.1    coxme_2.2-18.1      
[19] survival_3.5-7       bdsmatrix_1.3-6      survminer_0.4.9     
[22] ggpubr_0.6.0         extrafont_0.19       bootpredictlme4_0.1 
[25] merTools_0.6.1       marked_1.2.6         R2ucare_1.0.2       
[28] unmarked_1.3.2       devtools_2.4.5       usethis_2.2.2       
[31] effects_4.2-2        carData_3.0-5        broom.mixed_0.2.9.4 
[34] tidybayes_3.0.6      rptR_0.9.22          gt_0.9.0            
[37] patchwork_1.1.3      ggthemes_4.2.4       colorBlindness_0.1.9
[40] standardize_0.2.2    parameters_0.21.2    partR2_0.9.1        
[43] arm_1.13-1           lme4_1.1-35.1.9000   Matrix_1.6-1.1      
[46] MASS_7.3-60          gss_2.2-7            Rmisc_1.5.1         
[49] plyr_1.8.9           lattice_0.21-9       RColorBrewer_1.1-3  
[52] pbapply_1.7-2        BaSTA_1.9.5          readxl_1.4.3        
[55] lubridate_1.9.3      forcats_1.0.0        stringr_1.5.0       
[58] dplyr_1.1.3          purrr_1.0.2          readr_2.1.4         
[61] tidyr_1.3.0          tibble_3.2.1         ggplot2_3.4.3       
[64] tidyverse_2.0.0      RMark_3.0.0         

loaded via a namespace (and not attached):
  [1] estimability_1.4.1       msm_1.7                  coda_0.19-4             
  [4] knitr_1.44               multcomp_1.4-25          data.table_1.14.8       
  [7] useful_1.2.6             generics_0.1.3           snow_0.4-4              
 [10] leaflet_2.2.0            callr_3.7.3              terra_1.7-46            
 [13] cowplot_1.1.1            TH.data_1.1-2            proxy_0.4-27            
 [16] future_1.33.0            tzdb_0.4.0               webshot_0.5.5           
 [19] xml2_1.3.5               httpuv_1.6.11            xfun_0.40               
 [22] jquerylib_0.1.4          hms_1.1.3                ggdist_3.3.0            
 [25] satellite_1.0.4          evaluate_0.23            promises_1.2.1          
 [28] fansi_1.0.5              km.ci_0.5-6              DBI_1.1.3               
 [31] htmlwidgets_1.6.2        tensorA_0.36.2           stats4_4.2.2            
 [34] ellipsis_0.3.2           crosstalk_1.2.0          backports_1.4.1         
 [37] pbivnorm_0.6.0           insight_0.19.5           survey_4.2-1            
 [40] vctrs_0.6.4              abind_1.4-5              cachem_1.0.8            
 [43] withr_2.5.2              checkmate_2.2.0          emmeans_1.8.9           
 [46] prettyunits_1.2.0        mnormt_2.1.1             svglite_2.1.1           
 [49] crayon_1.5.2             leaflet.providers_1.13.0 labeling_0.4.3          
 [52] R2admb_0.7.16.3          pkgconfig_2.0.3          units_0.8-4             
 [55] nlme_3.1-163             pkgload_1.3.3            blme_1.0-5              
 [58] nnet_7.3-19              rlang_1.1.2              globals_0.16.2          
 [61] lifecycle_1.0.4          miniUI_0.1.1.1           extrafontdb_1.0         
 [64] cellranger_1.1.0         distributional_0.3.2     datawizard_0.9.0        
 [67] raster_3.6-23            KMsurv_0.1-5             zoo_1.8-12              
 [70] base64enc_0.1-3          processx_3.8.2           png_0.1-8               
 [73] KernSmooth_2.23-22       classInt_0.4-10          brew_1.0-8              
 [76] parallelly_1.36.0        rstatix_0.7.2            gridGraphics_0.5-1      
 [79] ggsignif_0.6.4           scales_1.2.1             memoise_2.0.1           
 [82] compiler_4.2.2           cli_3.6.1                urlchecker_1.0.1        
 [85] listenv_0.9.0            ps_1.7.5                 TMB_1.9.6               
 [88] tidyselect_1.2.0         stringi_1.7.12           mitools_2.4             
 [91] yaml_2.3.7               svUnit_1.0.6             survMisc_0.5.6          
 [94] sass_0.4.7               tools_4.2.2              timechange_0.2.0        
 [97] matrixcalc_1.0-6         uuid_1.1-1               rstudioapi_0.15.0       
[100] foreach_1.5.2            leafpop_0.1.0            snowfall_1.84-6.2       
[103] gridExtra_2.3            posterior_1.4.1          farver_2.1.1            
[106] digest_0.6.33            shiny_1.7.5              pracma_2.4.2            
[109] quadprog_1.5-8           Rcpp_1.0.11              car_3.1-2               
[112] broom_1.0.5              later_1.3.1              optimx_2023-8.13        
[115] sf_1.0-14                colorspace_2.1-0         fs_1.6.3                
[118] truncnorm_1.0-9          splines_4.2.2            expm_0.999-7            
[121] systemfonts_1.0.4        sessioninfo_1.2.2        xtable_1.8-4            
[124] jsonlite_1.8.7           nloptr_2.0.3             leafem_0.2.3            
[127] R6_2.5.1                 profvis_0.3.8            pillar_1.9.0            
[130] htmltools_0.5.6          mime_0.12                glue_1.6.2              
[133] fastmap_1.1.1            minqa_1.2.6              class_7.3-22            
[136] codetools_0.2-19         pkgbuild_1.4.2           mvtnorm_1.2-3           
[139] furrr_0.3.1              utf8_1.2.4               numDeriv_2016.8-1.1     
[142] arrayhelpers_1.1-0       curl_5.1.0               Rttf2pt1_1.3.12         
[145] CompQuadForm_1.4.3       rmarkdown_2.25           munsell_0.5.0           
[148] e1071_1.7-13             iterators_1.0.14         gtable_0.3.4            
[151] bayestestR_0.13.1